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Abstract 

Gain  scheduling,  the  traditional  method  of  providing  adaptive  control  to  a  nonlinear 
system,  has  long  been  an  ad  hoc  design  process.  Until  recently,  little  theoretical  guidance 
directed  this  practioners’  art.  For  this  reason,  a  systematic  study  of  this  design  process  and  its 
potential  for  optimization  has  never  been  accomplished.  Additionally,  the  nonlinearities  and 
the  large  seai'ch  space  involved  in  gain  scheduling  also  precluded  such  an  optimization  study. 
Traditionally,  the  gain  scheduling  process  has  been  some  variation  of  a  linear  interpolation 
between  discrete  design  points.  By  using  powerful  non-traditional  optimization  tools  such  as 
genetic  algorithms  there  aie  ways  of  improving  this  design  process. 

This  thesis  utilizes  the  power  of  genetic  algorithms  to  optimally  design  a  gain  sched¬ 
ule.  First,  a  design  methodology  is  validated  on  a  simple  pole  placement  problem,  then 
demonstrated  for  an  F-18  Super-maneuverable  Fighter.  From  this  experience,  a  general  gain 
scheduling  design  process  is  developed  and  presented. 


A  Gain  Scheduling  Optimization  Method 
Using  Genetic  Algorithms 


I.  Introduction 

1.1  Background 

Control  design  and  optimization  is  based  upon  many  assumptions  that  simplify  a  real 
world  control  problem  to  a  mathematically  manageable  model.  After  designing  a  controller 
for  a  simplified  model,  it  is  implemented  on  the  real  world  problem.  Consequently,  because 
of  the  simplifications  in  model  reduction,  there  are  difficulties  in  implementing  the  controller 
and  commanding  a  desired  response  throughout  the  entire  operating  envelope  of  the  system. 
One  method  of  overcoming  these  difficulties  is  gain  scheduling.  Gain  scheduling  is  a  method 
of  changing  the  controller  depending  upon  the  operating  condition.  There  are  four  basic  steps 
in  designing  a  gain  schedule  [30]: 

1.  Select  a  family  of  constant  operating  point  plants, 

2.  Design  a  controller  for  each  plant  in  family, 

3.  Design  a  method  of  scheduling  the  controllers  such  that  performance  is 
maintained  at  each  design  point,  and 

4.  Check  non-local  performance  of  the  scheduled  controller  by  simulation. 

Typically,  gain  scheduling  is  some  form  of  linear  interpolation  between  the  controllers  designed 
at  discrete  operating  points.  There  are  two  major  concerns  in  gain  scheduling.  The  primary 
concern  is  that  the  scheduled  controller  maintains  stability  throughout  the  operating  envelope 
[43,  46,  45,  44].  The  secondary  concern  is  that  acceptable  performance  is  also  maintained 
throughout  the  operating  envelope.  In  designing  a  gain  schedule,  there  are  many  design 
variables  the  control  designer  must  choose  [41]: 

•  The  variable  to  schedule  the  controller. 
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•  The  number  of  members  in  the  family  of  plants, 

•  The  members  in  the  family  of  constant  operating  point  plants, 

•  The  method  of  scheduling  the  controller, 

•  The  number  of  discrete  points  within  the  scheduling  variable  range  for  linear 
interpolation,  and 

•  The  distance  between  each  chosen  discrete  point  in  the  scheduling  variable 
range. 

In  numerically  optimizing  a  gain  schedule  design,  only  the  last  two  of  these  design  variables 
are  alterables.  Unfortunately,  traditional  calculus-based  optimization  methods  become  numer¬ 
ically  inaccurate  and  fail  to  find  an  optimal  solution  when  trying  to  simultaneously  optimize 
these  variables.  One  method  of  overcoming  this  difficulty  is  to  set  one  variable  and  optimize 
the  other,  however,  this  procedure  is  time  consuming  and  inefficient  when  there  are  means  of 
optimizing  these  variables  simultaneously.  Another  significant  difficulty  is  that  the  number 
of  the  design  variables  is  a  function  of  one  of  the  variables.  For  instance,  if  three  intervals  of 
the  scheduling  variable  are  selected  for  optimization,  then  there  are  not  as  many  variables  to 
optimize  as  if  eight  intervals  are  selected.  Additionally,  the  nonUnearities  of  the  system  and  the 
high  dimensionality  of  the  problem  makes  this  problem  very  difficult  for  traditional  calculus- 
based  deterministic  search  algorithms  [27,  32,  35].  Since,  genetic  algorithms  can  overcome 
these  difficulties  they  are  used  to  simultaneously  optimize  these  variables  [27,  32,  35]. 

Genetic  algorithms  (GAs)  have  found  global  optima  in  problems  with  both  nonlinearities 
and  high  dimensionality  where  calculus-based  methods  have  fallen  into  local  optimum,  such  as 
aircraft  structural  parameter  design  [9]  and  various  control  optimizations  [27, 32, 35, 38, 39]. 
For  these  reasons,  GAs  are  used  as  an  optimization  tool  in  this  investigation. 

1.2  Purpose 

The  purpose  of  this  research  effort  is  to: 

•  Demonstrate  that  gain  scheduled  controllers  can  be  optimized  using  GAs, 

•  Illustrate  the  feasibility  gain  scheduling  optimization  process  on  a  full  en¬ 
velope  aircraft  flight  control  example,  and 

•  Present  a  general  gain  schedule  design  method. 
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1.3  Scope  of  Research 

The  first  basic  step  in  gain  scheduling,  selection  of  operating  points,  is  a  procedure  that 
relies  on  common  sense  and  engineering  judgement,  consequentially,  it  is  not  reviewed  in  this 
research  effort.  Similarly,  the  second  basic  step,  controller  design,  is  independent  of  the  gain 
scheduling  design  process,  consequentially,  it  also  is  not  considered  in  this  research  work  [30]. 
Conversely,  the  third  step,  controller  scheduling,  can  be  optimized  with  respect  to  a  specified 
objective.  Therefore,  this  research  effort  is  limited  to  optimizing  this  procedure.  Finally,  the 
fourth  basic  step,  non-local  performance  validation,  merely  verified  the  global  performance 
of  the  designed  controller.  Additionally,  this  research  does  not  focus  on  improvements  in  the 
computation  efficiency  of  the  GA  optimization  algorithm  or  the  objective  function  calculation. 
This  effort  is  concerned  with  the  optimization  results  and  justifying  the  approach. 

1.4  Objectives 

The  specific  objectives  of  this  research  are; 

•  Find  a  readily  available,  easy  to  use  genetic  algorithm  platform  that  is 
compatible  with  a  common  computer-aided  control  design  tool, 

•  Formulate  an  appropriate  method  for  gain  schedule  optimization  while 
considering  alternatives, 

•  Validate  the  gain  scheduling  GA  optimization  process  on  a  simple,  repre¬ 
sentative  control  problem, 

•  Demonstrate  the  GA  optimization  process  on  a  full  envelope  flight  control 
example,  and 

•  Compile  the  lessons  learned  into  a  concise  design  optimization  procedure 
using  GAs. 

1.5  Preview 

The  next  two  chapters  provide  a  brief  summary  of  recent  published  work  concerning  gain 
scheduling  and  genetic  algorithms.  Chapter  II  summarizes  developments  in  the  theoretical 
aspects  of  gain  scheduling.  Chapter  III  presents  some  recent  GA  applications  in  the  field 
of  aircraft  controls  and  briefly  explains  how  GAs  work.  Additionally,  some  modifications 
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and  improvements  on  the  simple  GA  are  discussed,  and  finally  a  brief  GA  software  platform 
review  is  presented.  Chapter  IV  presents  a  simple  control  problem  to  validate  the  use  of  GAs 
in  gain  schedule  optimization.  Chapter  V  presents  information  for  an  F-18  gain  scheduled 
controller  designed  in  [4]  and  the  results  of  various  optimizations  of  this  gain  schedule. 
Chapter  VI  presents  a  general  gain  scheduling  optimization  design  process  developed  from 
the  experience  of  this  research  effort,  summarizes  the  results  of  this  research  work,  and 
provides  recommendations  for  future  research. 
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11.  Gain  Scheduling 


This  chapter  presents  a  review  of  gain  scheduling  and  how  it  has  been  achieved  in 
the  past.  By  understanding  the  methods  involved  in  scheduling  a  controller  throughout  an 
operating  envelope,  it  becomes  clearer  how  this  proven  technique  can  be  improved.  Section 
2.1  provides  an  explanation  of  why  gain  scheduling  is  used.  Next,  section  2.2  describes 
the  traditional  methods  of  gain  scheduling  which  are  improved  upon  in  this  research  effort. 
Concluding  this  chapter,  section  2.3  explains  how  a  gain  schedule  is  “optimized”  in  this 
investigation. 

2 .1  Rational  for  Gain  Scheduling 

Many  real  world  systems  are  nonlinear.  Also,  their  dynamics  may  change  as  a  function 
of  the  system’s  operating  condition.  There  are  two  basic  approaches  to  designing  a  controller 
for  these  systems.  The  preferred  way  is  to  design  a  controller  that  is  robust  enough  to  maintain 
stability  and  desired  performance  throughout  the  entire  operating  regime.  There  are  several 
techniques  for  designing  robust  controllers  including  Hoo  and  /^-synthesis  [14, 15].  This  can 
be  accomplished  for  some  problems  as  demonstrated  by  Shamma  and  Cloutier  in  designing 
a  missile  autopilot  [47].  Unfortunately,  this  is  not  always  possible.  The  second  approach  is 
to  change  a  linear  controller  as  the  plant  dynamics  change.  Hence,  the  controller  parameters 
(typically  gains)  are  scheduled  as  a  function  of  system’s  operating  condition.  Gain  scheduling 
has  proven  to  be  a  very  successful  method  of  implementing  a  global  controller  from  a  set  of 
linearized  controllers  [45]. 

There  are  some  significant  advantages  to  gain  scheduling  as  pointed  out  by  Rugh  [41]. 
The  main  advantage  is  that  linear  design  techniques  can  be  utilized.  When  designing  linear 
controllers  a  control  designer  has  a  wealth  of  computational  tools,  performance  measures, 
experience,  and  knowledge  to  draw  from  to  guide  the  design  process.  Another  advantage  is 
that  a  gain  scheduled  controller  has  the  potential  to  respond  rapidly  to  changing  operating 
conditions  whereas  more  modem  techniques  require  more  real  time  computation.  However, 
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there  are  some  difficulties  involved  in  the  design  process  [41].  The  major  difficulty  is  the 
selection  of  the  scheduling  procedure,  the  process  by  which  the  scheduling  variables  are 
changed.  Another  difficulty  is  in  the  selection  of  the  scheduling  variable.  Fortunately,  a 
couple  “rules  of  thumb”  have  emerged  to  help  overcome  this  second  difficulty.  First,  schedule 
on  a  variable  that  captures  the  nonlinearities  of  the  system,  and  second,  schedule  on  a  variable 
that  varies  slowly.  However  recent  work  by  Shamma  has  shown  that  these  rules  actually  have 
a  rigorous  mathematical  justification  [43, 46, 45, 44]. 

2.2  Gain  Scheduling  Process 

Following  is  a  typical  method  for  designing  a  gain  schedule  [44]. 

1.  Select  several  operating  points  covering  the  range  of  the  plant  dynamics. 

2.  Construct  a  linear  time  invariant  (LTI)  approximation  to  the  plant  and  design 
a  linear  controller. 

3.  Interpolate  the  controller  parameters  between  the  operating  points  in  order 
to  determine  control  parameters  for  operating  points  between  the  selected 
design  operating  points. 

To  interpolate  the  controller  parameters  simple  curve  fitting  techniques  are  used  [41,  43]. 
However,  there  are  no  guarantees  of  stability  or  robustness  without  numerous  simulations 
or  theoretical  analysis.  An  idealized  gain  schedule  is  defined  as  a  gain  schedule  that  has  a 
controller  designed  at  every  operating  point  within  the  operating  envelope  [41  ].  Obviously,  for 
any  large  number  of  operating  points  the  idealized  gain  schedule  is  impractical.  Nonetheless, 
approximations  to  this  idealized  gain  schedule  can  be  done  where  the  resulting  controller  is 
selected  by  a  table-look  up  method. 

Fortunately,  several  methods  and  theorems  have  been  published  to  analyze  the  system 
for  global  stability  [29, 43, 46, 45, 44].  Unfortunately,  these  methods  are  based  on  scheduling 
simple  controllers  such  as  proportional-integral-derivative  PID  controllers  and  full  state  feed 
back  gain  matrices.  These  methods  do  not  directly  address  the  complexities  of  modern 
controller  designs  such  as  H2,  and  n  synthesis. 
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Recent  works  have  been  focusing  on  scheduling  these  modem  controllers  [6, 5, 18, 37]. 
These  works  use  a  method  of  linear  fractional  transforms  (LFTs)  to  approach  to  problem  of 
gain  scheduling.  The  scheduling  variable  is  included  as  an  exogenous  signal  to  be  tracked  or  as 
a  disturbance  to  be  rejected.  A  bounded  controller  space  is  solved  via  linear  matrix  inequality 
(LMI)  equations.  Using  this  method,  a  controller  for  a  linear  parameter  varying  (LPV)  plant 
is  linear  parameter  varying  itself.  The  resulting  controller  is  thereby  self  scheduled. 

2.3  Gain  Scheduling  Design 

The  traditional  approach  to  gain  scheduling  is  selected  because  of  its  wide  use  and 
acceptance.  The  goal  of  the  design  is  to  directly  optimize  the  controller  with  respect  to  an 
objective  function  over  the  operating  range  of  the  system.  This  is  different  than  optimizing  a 
controller  at  a  specific  operating  point  and  then  finding  a  linear  fit  of  the  controllers  over  the 
operating  range.  Ideally  the  control  designer  would  like  to  design  a  controller  that  responds 
exactly  the  same  regardless  of  the  operating  condition.  A  gain  scheduled  controller  of  this 
type  would  be  an  optimal  gain  scheduled  controller  defined  as  follows: 

Definition  1  (Optimal  Gain  Scheduled  Controller)  An  optimal  gain  scheduled  controller 
provides  a  uniform  response  throughout  the  operating  envelope. 

The  optimality  of  a  gain  schedule  can  be  measured  by  comparing  a  desired  closed  loop  response 
to  an  actual  response  at  a  specific  operating  point.  The  desired  response  can  be  defined  in 
either  the  time  or  frequency  domain.  The  objective  of  the  design  optimization  is  to  minimize 
the  deviation  between  the  actual  response  and  the  desired  response  at  various  points  in  the 
operating  envelope.  In  other  words,  the  gain  scheduling  error  is  minimized.  Following  are 
definitions  to  clarify  the  difference  between  tracking  error  and  the  gain  scheduled  error  that  is 
minimized  in  this  investigation. 

Definition  2  (Controller  Tracking  Error)  The  difference  between  the  commanded  input  and 
the  closed  loop  output  response  at  one  specific  operating  point. 
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Definition  3  (Gain  Scheduling  Error)  The  difference  between  the  closed  loop  response  of 
the  central  controller  and  a  non  central  controller. 

Definition  4  (Central  Controller)  The  controller  selected  at  a  specific  operating  point  as 
having  the  desired  response.  The  closed  loop  response  of  this  controller  is  used  in  calculating 
the  gain  scheduling  error. 

Therefore,  the  specific  goal  of  the  gain  schedule  design  optimization  is  to  minimize  the 
gain  scheduling  error  in  the  closed  loop  response  throughout  the  operating  envelope.  The 
central  controller  can  be  a  set  of  desired  eigenvalues,  a  set  of  dominant  closed  loop  poles,  or 
a  controller  designed  at  a  specific  operating  point  [4].  The  example  optimization  discussed 
in  chapter  IV  selects  a  set  of  dominant  closed  loop  poles  to  be  the  central  controller,  and  in 
chapter  V  the  central  controller  is  a  controller  defined  at  a  specific  operating  point.  Each 
optimization  performed  demonstrates  how  the  gain  scheduling  error  can  be  minimized. 

To  clarify  the  difference  the  scheduling  variable  and  the  scheduled  variable,  they  are 
defined  as  follows: 

Definition  5  (Scheduling  Variable)  The  variable(s)  that  is  representative  of  the  changes  in 
the  plant  dynamics  and  is  a  measurable  signal. 

Definition  6  (Scheduled  Variable)  The  variable(s)  that  are  changed  as  a  function  of  the 
scheduling  variable,  thereby  changing  the  controller  as  a  function  of  the  operating  point. 

2 .4  Optimization  Approach 

The  first  step  is  to  discretize  the  scheduling  variable  range  by  breaking  it  into  N 
intervals.  At  first  the  number  of  intervals  are  constrained  to  a  arbitrary  number.  Additionally, 
the  size  of  the  intervals  are  also  constrained,  for  simplicity  all  constrained  interval  sizes  will 
be  equal.  Next,  through  some  preliminary  analysis  the  scheduled  variables  are  bounded 
by  their  minimum  and  maximum  values  over  the  scheduling  variable  range.  After  defining 
an  appropriate  objective  function,  the  gain  scheduling  error  is  minimized  by  scheduling  the 
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controller  parameters  as  piecewise  constant  or  piecewise  linear  functions  of  the  scheduling 
variable.  For  each  interval,  the  scheduled  variables  are  allowed  to  vary  independently  between 
their  minimum  and  maximum  bounds.  Further  optimization  unconstrain  the  size  of  the 
intervals  and  the  number  of  intervals  used  to  discretize  the  scheduling  variable. 
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III.  Genetic  Algorithms 


Genetic  Algorithms  (GAs)  are  an  optimization  method  based  on  Darwin’s  theory  of 
“Survival  of  the  Fittest.”.  GAs  were  first  developed  by  Holland  in  1975  [25].  Since  then 
they  have  been  extensively  expanded  by  DeJong  [1 1, 12],  Goldberg  [20, 21, 22],  Grefenstette 
[23]  and  many  others.  GAs  differ  from  classical  methods  of  optimization  in  many  ways.  The 
next  section  describes  these  differences.  The  following  two  sections  briefly  explain  how  and 
why  GAs  work.  Some  improvements  on  these  simple  GAs  are  presented  in  section  3.4.  A 
review  of  GA  software  is  presented  in  section  3.5.  Section  3.6  concludes  this  chapter  with  an 
explanation  of  the  implementation  of  the  GA. 

3.1  Classical  Optimization 

Classical  methods  of  optimization  fall  into  three  categories:  calculus-based,  enumera- 
tive,  and  random.  Consider  a  general  optimization  problem  statement  [50]: 


Minimize:  F{X)  the  objective  function  (3.1) 

Subject  to: 

gj(X.)  <0  j  =  l,m  inequality  constraints  (3.2) 

hkCX.)  =  0  k=l,l  equality  constraints  (3.3) 

<  Xj  <  X“  i  =  l,n  side  constraints  (3.4) 

X2 

where  X  =  <  X3  >  the  design  variables 


I  J 
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Calculus  based  methods  optimize  this  given  problem  by  satisfying  three  Kuhn-Tucker  condi¬ 
tions  [50].  The  Kuhn-Tucker  conditions  are  (X*  is  the  optimal  design  vector): 


1.  X*  is  feasible 


2.  =  0  i  =  l,m 

m 

3.  Vf(X*)+5:A^Vft(X-) 

i=i 

A,>0 


where  VF(X)  = 


A,>0 

i 

+  Y.^k+mXhk{X*)  =  0 

k=l 

A  jt+m  unrestricted 
6F{X)/SXi 
6F{X)/6X2 
6F{X)/6X3 


(3.5) 

(3.6) 

(3.7) 


(3.8) 


6F{X)/8Xn 


Both  the  power  and  the  drawback  of  calculus-based  methods  becomes  apparent  in  the  third 
Kuhn-Tucker  condition,  Eq  (3.7).  The  gradient  of  both  the  objective  function  and  the  con¬ 
straints  must  be  evaluated  to  determine  the  direction  the  search  wiU  proceed.  For  a  function 
where  these  gradients  do  not  have  analytical  solutions,  they  must  be  evaluated  by  a  finite- 
difference  method  [50].  Additionally,  for  the  calculus-based  algorithm  to  determine  whether 
the  optimum  that  it  has  found  is  a  maximum  or  minimum,  it  must  evaluate  the  Hessian  matrix 
of  second  derivatives  (Eq  (3.9)). 


S^F(X) 

6^F(X) 

S^FiX) 

SX'f 

5X1SX2 

6Xi6X„ 

6^FiX) 

S^F{X) 

H  = 

SX26X1 

^xT 

SX2SX„ 

5^F{X) 

52F(X) 

S^F{X) 

L  SXnSXi 

SX„SX2 

sxi  J 

Even  though  calculus-based  methods  can  be  very  efficient,  they  “break  down”  when  the 
objective  function  or  its  gradient  is  not  continuous  or  “well  behaved”.  Moreover,  the  solution 
is  dependent  upon  the  initial  conditions  given  to  the  algorithm.  For  a  given  function,  a 
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calculus-based  algorithm  is  guaranteed  to  find  the  nearest  local  extremum  if  it  is  exists.  The 
solution  is  a  local  minimum  if  the  Hessian  matrix  is  positive  definite  at  the  solution  point  [50]. 
Unfortunately,  the  solution  found  may  still  not  be  the  global  minimum. 

Enumerative  techniques  were  developed  to  guarantee  that  the  global  minimum  was 
found  by  exhaustively  searching  all  points  in  the  design  space.  Unfortunately,  as  the  number 
of  design  variables  increases  the  number  of  possible  solutions  increases  exponentially.  This 
is  the  “curse  of  dimensionality.”  Obviously  for  problems  with  a  large  number  of  design 
variables,  the  enumerative  technique  is  inefficient. 

Random  methods  were  designed  to  increase  the  efficiency  of  the  enumerative  methods 
by  randomly  evaluating  points  in  the  search  space.  Since  only  random  points  in  the  solution 
space  are  evaluated  there  is  no  guarantee  of  finding  the  global  optimum,  therefore  the  algorithm 
continues  until  a  large  percentage  of  the  search  space  has  been  evaluated.  Again,  as  the  number 
of  design  variables  increases  the  number  of  possible  solutions  increases  at  a  greater  rate.  Since 
random  methods  lack  efficiency,  they  are  not  expected  to  obtain  good  results  [27  ] .  G As  provide 
an  approach  that  overcomes  these  difficulties. 

GAs  differ  from  traditional  optimization  techniques  in  four  basic  ways:  1)  they  code 
the  parameter  set  instead  of  the  parameters  themselves,  2)  they  search  a  population  of  points 
instead  of  a  single  point  in  the  solution  space,  3)  they  do  not  require  an  initial  guess  of  the 
solution,  and  4)  they  use  probabilistic  transition  rules  instead  of  deterministic  transition  rules. 
GAs  are  a  zero  order  optimization  method  that  does  not  depend  on  the  behavior  of  the  objective 
function.  Consequently,  they  have  provided  robust  optimization  over  a  wide  range  of  problems 
[3, 27, 32].  Some  problems  that  have  been  optimized  using  GAs  are  shipping  scheduling  [34], 
dynamic  control  [32],  control  optimization  [35,  27],  floor  plan  area  optimization  [40],  and 
investment  market  timing  strategies  [7].  Traditional  and  genetic  optimization  techniques  were 
compared  for  efficiency  and  accuracy  in  an  aircraft  parameter  design  optimization  [9].  In 
the  study  the  authors  performed  an  extensive  comparison  of  calculus,  enumerative,  random, 
genetic,  and  simulated  annealing  algorithms.  GAs  and  Simulated  Annealing[26]  were  the 
only  two  methods  to  find  the  global  optimum  out  of  the  fifteen  algorithms  tried. 
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Additionally  advantages  of  GAs  are:  1)  they  do  not  require  the  derivative  of  the  function 
to  exist,  2)  they  can  handle  objective  functions  with  penalty  functions  without  difficulty,  unlike 
a  calculus-based  method  that  must  transform  the  objective  function  into  a  more  manageable 
form  [27],  3)  as  the  dimension  of  the  problem  grows  the  computation  time  for  the  GA  grows 
in  a  linear  fashion,  whereas  calculus-based  methods  grow  at  least  quadratically  [34,  35],  4) 
they  are  not  dependent  on  the  shape  of  the  solution  space  for  finding  an  optimum,  and  5)  the 
optimum  found  by  the  GA  is  more  likely  to  be  the  global  optimum  value  for  the  objective 
function  since  the  solution  is  not  dependent  upon  the  initial  condition  given  to  the  algorithm. 

However,  there  are  a  couple  of  disadvantages  to  GAs.  First,  GAs  are  a  naturally  parallel 
algorithm  that  are  typically  run  on  a  serial  machine  [32].  This  observation  was  also  made  by 
Goldberg  [20]: 

“In  a  world  where  serial  algorithms  are  usually  made  parallel  through  countless 
tricks  and  contortions,  it  is  no  small  irony  that  genetic  algorithms  are  made  serial 
through  equally  unnatural  tricks  and  turns.” 

The  inefficiency  of  serial  computation  significantly  increases  the  computation  time  of  the 
GA.  Additionally,  parallel  GAs  require  fewer  function  evaluations  than  the  best  alternative 
algorithms  [32].  Another  disadvantage  of  GAs  is  that  they  are  robust  methods  for  searching 
the  solution  space  to  find  the  area  in  which  the  global  optimum  occurs,  but  they  are  inefficient 
at  obtaining  a  solution  where  high  precision  is  desired.  The  following  paragraph  reviews  some 
literature  that  compares  GAs  and  traditional  calculus-based  algorithms. 

A  GA  was  compared  to  Powell’s  conjugate  search  direction  for  two  control  system 
optimization  problems  [27].  The  first  problem  was  a  linear  quadratic  regulator  problem  for 
a  lateral  autopilot  which  was  solved  via  the  algebraic  Riccati  equation  [19]  for  the  exact 
answer.  Both  the  GA  and  Powell’s  method  found  the  exact  solution.  A  significant  difference 
was  that  Powell’s  method  required  an  order  of  magnitude  more  function  evaluations.  A 
second  comparison  was  performed  on  a  wind  shear  problem  where  the  objective  was  to  find 
a  controller  that  minimized  the  deviation  in  the  velocity  and  flight  path  of  an  aircraft  due  to  a 
wind  burst.  The  GA  was  able  to  find  the  exact  solution,  but  Powell’s  method  was  dependent 
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upon  the  initial  conditions  and  unable  to  find  the  global  optimum.  Similar  experiments  were 
performed  by  Michalewicz  [34],  who  drew  the  following  conclusions  on  the  effectiveness  of 
GAs  applied  to  optimal  control:  1)  the  objective  function  need  not  be  continuous  over  the 
search  space  for  GAs  to  find  an  optimum;  2)  GAs  give  intermediate  information  while  the 
problem  is  being  solved  so  that  when  a  desired  degree  of  accuracy  is  reached  computation 
can  be  halted,  whereas  other  packaged  optimization  methods  do  not  return  a  result  until  the 
optimization  is  complete;  and  3)  GAs  grow  linearly  with  the  size  of  the  problem.  GAs  are  not 
the  best  optimization  method  for  all  problems,  but  for  problems  with  high  dimensionality  and 
nonlinearities,  GAs  do  have  an  advantage. 

3.2  GAs:  The  Inner  Workings 

GAs  explore  and  exploit  various  solutions  by  manipulating  building  blocks  called 
schemata.  GAs  explore  different  possible  solutions,  and  exploit  the  best  solutions  through 
the  implementation  of  three  genetic  OTpemots-rnutation,  crossover,  and  reproduction.  GAs 
evolve  a  population  of  solutions  from  generation  to  generation  through  these  genetic  operators 
to  an  optimal  solution. 

GAs  can  either  minimize  or  maximize  a  given  objective  function.  Any  minimization 
problem  can  be  transformed  into  an  equivalent  maximization  problem  through  the  mapping 
depicted  in  Eq  (3.10). 

Qmax  —  fmin  (3.10) 

Therefore,  for  simplicity  all  problems  in  this  chapter  will  be  considered  maximization  prob¬ 
lems. 

GAs  encode  the  parameter  range  of  the  design  variable  and  map  it  to  a  binary  string 
of  a  specified  length.  For  example,  a  binary  string  of  length  seven  has  a  minimum  value, 
0000000,  which  is  mapped  to  the  minimum  value  of  the  design  variable,  x^in-  The  maximum 
value  of  the  string,  1111111,  is  mapped  to  the  maximum  value  of  the  design  variable,  x^ax- 
Between  these  two  extremes  there  is  a  linear  mapping  of  the  string  values  to  the  parameter 
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values.  For  a  function  of  more  than  one  variable,  the  strings  representing  each  variable  are 
concatenated  together  to  form  one  chromosome.  For  example,  a  function  with  two  design 
variables,  with  each  having  a  string  length  of  seven  would  form  a  chromosome  of  length 
fourteen,  00000001111111. 

The  following  two  sections  will  explain  how  the  GA  explores  and  exploits  solutions 
through  genetic  operators. 

3.2.1  Exploration.  The  main  operator  responsible  for  exploration  of  the  search 
space  is  mutation.  The  mutation  operator  is  controlled  by  a  user  specified  probability,  Pm-  This 
probability  is  the  chance  that  each  bit  in  the  chromosome  string  is  mutated  to  its  complement 
value  of  0  or  1.  For  each  bit  in  the  population,  a  random  number  between  zero  and  one  is 
generated.  If  the  random  number  is  less  than  pm ,  that  bit  is  mutated.  For  example,  suppose  that 
a  chromosome  has  the  form  shown  in  Fig.  3.1(a).  If  the  fourth  bit  is  chosen  for  mutation,  the 
chromosome  is  changed  to  the  new  chromosome  shown  in  Fig.  3.1(b).  This  new  chromosome 
represents  a  new  possible  solution. 


1011010110 

(a) 

1010010110 

(b) 

Figure  3.1  Mutation 

The  frequency  of  mutation  for  a  specified  probability,  Pm,  is  dependent  on  the  string 
length,  m,  of  the  chromosome.  Eq  (3.11)  describes  the  expected  number  of  mutations  per 
generation. 

E[#  of  mutations]  =  ny.my.pm  (3.11) 

where  n  =  the  population  size,  m  =  the  string  length,  and  Pm  =  the  user  specified 
mutation  probability. 

For  a  string  length  m  =  100  and  pm  =  0.01,  the  expected  number  of  mutations  per  string  is 
1.  For  a  string  length  m  —  50,  then  expected  number  of  mutations  is  only  1  per  every  two 
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01110110110 
10101101100 
(a)  Parents 
01110101100 
10101110110 
(b)  Offspring 

Figure  3.2  Crossover 

strings  in  the  population.  Consequently,  for  a  population  size  of  n  =  20,  a  string  length  of 
m  —  50,  and  Pm  =  0.01,  the  expected  number  of  mutations  per  generation  is  10.  Typically, 
values  of  Pm  range  from  0.01  to  0.001  [23].  For  values  of  pm  >  0.05,  the  GA  approaches  a 
random  search  method. 

3.2.2  Exploitation.  Two  genetic  operators  are  responsible  for  exploitation  of 
good  solutions,  crossover  and  reproduction.  The  following  two  sections  wUl  describe  these 
operators  in  more  detail. 

3. 2. 2.1  Crossover.  The  crossover  operator  acts  on  two  solution  strings  si¬ 
multaneously  by  swapping  information  included  in  each  string.  The  probability  that  crossover 
occurs.  Pc,  is  set  by  the  user.  Once  two  strings  are  chosen  for  crossover,  a  position  between 
the  bits  is  randomly  chosen  with  a  uniform  distribution.  Suppose  that  the  two  strings  shown 
in  Fig.  3.2(a)  are  chosen  for  crossover,  and  that  the  position  between  the  fifth  and  sixth  bit 
is  randomly  selected  as  the  crossover  point.  At  the  crossover  point  the  strings  are  cut.  The 
partial  strings  are  swapped,  thereby  forming  the  two  new  strings  shown  in  Fig.  3.2(b). 

Crossover,  unlike  mutation,  is  independent  of  the  string  length,  but  is  still  a  function  of 
the  population  size.  Eq  (3.12)  describes  the  expected  number  of  crossovers  per  generation. 

E[#  of  strings  selected  for  crossover]  =  n  x  pc  (3.12) 

For  a  population  size  of  n  =  50,  and  Pc  =  0.6,  the  expected  number  of  strings  selected  for 
crossover  would  be  30.  Typically,  values  for  pc  range  from  0.6  to  0.9  [23]. 
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Table  3.1  Reproduction  Example 


string 

fitness  value 

fiftot 

01101001 

28.6 

10010010 

16.4 

01011011 

23.7 

0.30 

10110110 

9.8 

0.13 

78.5 

1.00 

3 .2 .2 .2  Reproduction.  The  reproduction  operator  acts  on  the  population  only 
after  both  the  mutation  and  crossover  operators  have  been  implemented.  The  reproduction 
operator  reproduces  each  solution  with  a  probability  proportional  to  the  fitness  value  of  each 
string.  The  strings  with  above  average  fitness  values  are  reproduced  at  higher  rates  than  strings 
with  below  average  fitness  values.  Consequently,  poor  solutions  die  out  of  the  population.  By 
not  eliminating  the  worst  solutions  immediately,  the  population  remains  diverse  and  therefore 
does  not  converge  prematurely  to  a  local  optimum. 

Consider  a  population  of  size  n  =  4  shown  in  Table  3.1.  Each  string  is  evaluated 
using  an  objective  function  to  find  the  respective  fitness  of  each  string.  The  total  sum  of  all 
the  fitness  values  is  78.5.  Each  string  has  a  chance  of  being  selected  for  the  next  generation 
proportional  to  its  fitness  value,  /,  divided  by  the  total  fitness  value  of  the  population,  ftot- 
Therefore,  the  probability  that  a  string  is  selected  for  the  next  generation  is: 

/ 

^(selection  for  next  generation)  =  (3.13) 

Jtot 

The  sum  of  these  probabilities  totals  100  per  cent.  Reproduction  takes  place  by  giving  each 
string  a  space  on  a  roulette  wheel  proportional  to  its  relative  fitness,  f/ ftot  (See  Fig.  3.3). 
The  roulette  wheel  is  spun  four  times  to  select  the  individuals  for  the  next  generation. 

This  concludes  the  basics  of  a  GA.  The  next  section  will  explain  how  these  genetic 
operators  work  together  to  find  an  optimal  solution  to  a  given  objective  function. 
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30% 


Figure  3.3  Reproduction  Roulette  Wheel 


33  Schemata 

Schemata  ^  are  a  representation  of  a  set  of  binary  strings  used  in  the  theoretical  structure 
of  GAs.  To  understand  schemata,  a  “don’t  care”  symbol  is  introduced.  A  schema  of  length 
eight,  5o  =  {**101011},  would  represent  any  string  of  the  same  length  that  had  Is  and  Os  in 
the  same  positions  as  the  schema  itself.  A  schema  with  r  “don’t  care”  symbols  has  2’’  strings 
matching  that  schema.  Every  schema  matches  2’’  strings,  where  r  is  the  number  of  “don’t 
care”  symbols  in  the  schema.  Thus,  the  following  strings  would  match  schema  So. 

(00101011),  (01101011),  (10101011),  (11101011) 

Conversely,  each  string  of  length  m  is  matched  by  2’"  schemata.  For  example,  the  string 
(101 10)  is  matched  by  the  following  2®  schema: 

(10110) 

(*0110) 

(1*110) 

(**110) 

(10*10) 


^Note:  Schemata  is  the  plural  of  schema. 
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There  are  two  main  characteristics  used  to  define  schema,  order  and  defining  length. 
Order,  denoted  by  o(S),  is  the  number  of  Is  and  Os  specified  in  the  schema.  It  defines  the 
specialty  of  the  schema.  A  high  order  schema  is  more  specific  than  a  low  order  schema.  A 
schema  of  high  order  is  more  likely  to  be  disrupted  by  mutation  than  a  low  order  schema.  The 
order  of  two  schemata  are  shown  in  Fig.  3.4. 

The  defining  length,  denoted  by  S{S),  of  a  schema  is  the  distance  between  the  first 
and  last  fixed  positions.  The  defining  length  quantifies  the  compactness  of  a  schema.  The 
minimum  defining  length  is  zero  for  a  single  specified  position,  and  the  maximum  defining 
length  for  a  string  of  length  m  is  (m  - 1).  Defining  length  is  a  determining  factor  in  a  schema 
surviving  crossover,  with  compact  schemata  more  likely  to  survive.  The  defining  length  of 
two  schemata  are  also  shown  in  Fig.  3.4. 

5o={***ll***}  o(5o)  =  2  ,5(5o)  =  5-4  =  1 

5i={*10**l*l}  o(5i)=4  ^(5i)  =  8-2  =  6 

Figure  3.4  Order  and  Defining  Length 


33.1  Survival  of  the  Fittest.  For  a  population  of  strings  with  length  m,  there  are 
a  total  of  3™  possible  schema.  With  a  population  of  size  n,  there  will  only  be  2™  to  n  x  2”^ 
schemata  represented,  with  some  more  fit  than  others.  The  fitness  of  a  schema  is  determined 
by  the  average  fitness  of  the  strings  matching  schema  S.  Let  N(S,t)  represent  the  number  of 
strings  in  generation  t  that  match  schema  S.  Then  the  following  equation  calculates  the  fitness 
of  a  schema  S  in  generation  t. 

N(S,t)  r 

e.al(S,t)=  E  (3.14) 

From  this  relation  the  expected  number  of  strings  in  the  next  generation,  t  +  1,  can  be 
determined  for  a  population  with  given  average  fitness  (favg  =  ftotin). 

E[N{S,  t  +  1)]  =  N{S,  (3.15) 

Javg 
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A  schema  with  an  above  average  fitness  will  receive  an  increasing  number  of  strings  in  the 
subsequent  generations  as  shown  in  Eq  (3.15).  This  increase  is  exponential  for  a  schema  with 
a  consistently  above  average  fitness  [25, 34]. 

However,  the  genetic  operators,  crossover  and  mutation,  can  destroy  a  schema  and 
therefore  reduce  this  growth  rate  of  the  schema.  Consequently  the  probability  of  crossover 
destroying  a  schema  is; 

Pd  = - 7  X  Pc  (3.16) 

m  —  1 

Therefore  the  probability  of  survival  is: 

Ps  =  l-  X  p,  (3.17) 

However,  there  is  a  chance  that  even  if  crossover  occurs  within  the  defining  length  of  a  schema, 
the  schema  will  survive.  For  an  example  of  this  see  Fig.  3.5.  Therefore,  the  probability  of 
survival  is  slightly  larger: 

P.  >  1  -  X  Pe  (3.18) 

5'o={***ll***} 

011111011 
110111100 
(a)  Parents 
011111100 
110111011 
(b)  Offspring 

Figure  3.5  Schema  Surviving  Crossover 

The  probability  that  a  schema  survives  mutation  depends  on  its  specialty: 

Ps  =  (1  -  Pdl-praX  o{S)  for  Pm  <  1  (3.19) 
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Combining  the  effects  of  crossover  and  mutation,  Eq  (3.15)  becomes: 

£(iV(S,i  +  1)]  >  [l  -  Po-^  -  P„o(S)]  (3.20) 

Javg  Tn  1 

To  demonstrate  how  this  growth  works  consider  the  following  schema: 

g_|**U^Q*********!i:*!t;=t:**!t:| 

m  =  21 

6{S)  =  5-3  =  2 
o(5)  =  3 
N{S,t)  =  A 
eval{S,  t)  =  20 

favg  =  12 

Pc  =  0.6 
Pm  =  0.001 

£;[iV(5,f  +  1)]  >  4.0  X  1.667  x  0.938  w  6 

For  a  below  average  schema  the  number  of  strings  in  subsequent  generations  would 
decrease  in  a  similar  manner.  The  following  theorem  describes  the  concept  given  quantitatively 
in  Eq  (3.20).  This  theorem  is  known  as  the  fundamental  theorem  of  GAs  [34]. 

Theorem  1  (Schema  Theorem)  Short,  low-order,  above  average  schemata  receive  exponen¬ 
tially  increasing  trials  in  subsequent  generations  of  a  genetic  algorithm. 

3.4  Improvements  on  the  Simple  GA 

Several  modifications  have  been  proposed  to  improve  the  performance  of  the  simple 
GA  (sGA).  One  of  the  simplest  ways  of  improving  the  precision  of  the  GA  is  to  use  a  floating 
point  representation  of  the  design  variables  instead  of  mapping  them  to  binary  strings.  Two 
different  research  efforts  have  found  that  the  use  of  floating  point  genes  has  increased  the 
speed,  precision  and  the  accuracy  of  the  GA  [10, 34]. 
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There  are  two  main  qualities  of  the  GA  that  have  proven  it  to  be  a  robust  optimization 
method  for  a  wide  range  of  practical  problems:  its  power  to  explore  new  solutions,  and  the 
ability  to  exploit  the  best  solutions  it  has  found.  The  user  has  control  over  these  parameters  by 
selecting  the  crossover  probability,  Pc,  and  the  mutation  probability,  Pm-  For  improved  results 
the  user  might  want  to  vary  these  parameters  over  the  computation  time  of  the  algorithm.  The 
following  two  sections  discuss  some  of  the  ways  that  crossover  and  mutation  operators  may 
be  varied  or  adaptively  changed.  The  third  section  3.4.3  discusses  the  advantages  of  using  a 
messy  GA  (mGA)  instead  of  a  sGA. 

3.4.1  Crossover  Operator .  One  simple  way  of  changing  the  crossover  operator  is 
to  vary  the  number  of  points  at  which  crossover  can  occur.  By  allowing  multi-point  crossover 
to  occur,  some  schemata  with  long  defining  lengths  can  be  preserved  which  would  otherwise 
be  destroyed  with  single  point  crossover.  For  example,  two  point  crossover  would  randomly 
select  two  points  along  the  length  of  the  chromosome  for  crossover  (See  Fig.  3.6).  Then  the 
portion  of  the  string  between  the  two  crossover  points  are  swapped  between  the  two  parents 
to  produce  two  new  offspring. 


OlllOllOllOlllOlll 
lOlOllOllOOIOllOlO 
(a)  Parents 
OlllOIOllOOIllOlll 
lOlOlllOllOIOllOlO 
(b)  Offspring 

Figure  3.6  Two  Point  Crossover 

Further  generalization  of  this  concept  can  lead  to  the  idea  of  uniform  crossover.  Uniform 
crossover  decides  with  probability,  Pc,  which  bit  positions  of  the  first  parent  will  be  exchanged 
with  the  second  parent.  This  is  similar  to  the  mutation  operator  in  the  sense  that  each  bit  has  a 
chance  to  be  crossed  with  another  bit  from  the  second  parent.  For  an  example  of  this  type  of 
crossover,  see  Fig.  3.7  {pm  =  0.5). 
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1001100110 
0110101011 
(a)  Parents 
1100100110 
0011101011 
(b)  Offspring 

Figure  3.7  Uniform  Crossover 

Another  innovative  way  of  changing  the  crossover  operator  is  to  change  the  probability 
that  crossover  occurs  as  a  function  of  the  average  and  best  fitness  of  the  population  [48].  This 
operator  is  shown  in  Eq  (3.21). 


Pc  =  hi  f max  -  f)/ifn 
Pc  =  h 


favg)  for/'  >  favg 

for/'  <  favg 


(3.21) 


where  /'  is  the  largest  fitness  value  of  the  two  parents  selected  for  crossover, 
fmax  is  the  largest  fitness  value  of  the  population,  and  favg  is  the  average  fitness 
value  of  the  population. 

The  parameters  h  and  k2  are  scaling  parameters  that  can  be  chosen  arbitrarily.  The  authors 
of  [48]  chose  h  =  1-0  and  k2  =  1.0.  With  this  crossover  operator  each  individual  in  the 
population  has  a  different  crossover  probability  depending  on  its  fitness  value.  Individuals  with 
above  average  fitness  values  have  a  higher  crossover  probability  than  those  with  below  average 
fitness.  The  probability  of  crossover  is  0.0  for  individuals  with  fitness  values  equal  to  fmax- 
Davis  [10]  presented  another  method  of  adaptive  crossover  where  the  crossover  probability  is 
dependent  upon  the  fitness  of  the  offspring  produced.  The  greater  the  offspring’s  fitness,  the 
greater  the  probability  that  crossover  will  occur  at  that  bit  location. 

3.4.2  Mutation  Operator .  The  choice  of  Pm  is  critical  in  the  performance  of  the  G A 
[11].  Two  ways  of  adaptively  changing  the  mutation  operator  are  presented.  Both  have  shown 
to  provide  increased  performance  in  the  GA  in  a  sampling  of  test  problems  [34, 48].  The  first 
method  changes  Pm  as  a  function  of  the  current  generation  number  and  the  total  number  of 
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generations.  This  mutation  operator  was  designed  as  a  floating  point  mutation  operator,  but 
it  can  be  transformed  into  a  binary  operator  with  an  appropriate  mapping  [34].  This  mutation 
operator  is  shown  in  Eq  (3.22). 


(3.22) 

where  T  is  the  maximum  number  of  generations,  t  is  the  current  generation 
number,  r  is  a  random  number  G  [0,1],  6  is  a  system  parameter  determining  the 
degree  of  dependence  on  iteration  number  (a  value  of  5  is  used  for  experimental 
results),  and  y  is  the  value  of  the  gene  being  mutated. 

The  effect  of  this  operator  is  that  when  the  GA  begins,  there  is  a  high  probability  of  mutation. 
This  maintains  the  diversity  of  the  population  searching  for  new  solutions.  As  the  GA  nears 
completion,  t  T,  mutation  is  decreased  allowing  crossover  to  become  dominant  and 
converge  to  an  optimal  solution. 

Another  published  adaptive  mutation  operator  changes  the  mutation  probability  depend¬ 
ing  on  an  individual’s  fitness  value  instead  of  the  generation  number  [48].  Eq  (3.23)  defines 
the  proposed  mutation  operator. 

Pm  ~  f}/{fmax  favg)  for  f  ^  favg 

Pm  =  k4  for  /  <  favg  (3.23) 


where  /  is  the  fitness  of  the  individual,  fmax  is  the  maximum  fitness  in  the 
population,  and  favg  is  the  average  fitness  of  the  population. 

The  values  of  ks  and  k^  can  be  chosen  by  the  user;  the  authors  of  [48]  chose  values  of  0.5  for 
both.  Individuals  with  above  average  fitness  have  lower  mutation  probabilities  than  individuals 
with  below  average  fitness.  The  mutation  probability  is  0.0  for  individuals  with  fitness  equal 
to  fmax-  Individuals  with  below  average  fitness  are  totally  disrupted.  The  combined  effect  of 
this  mutation  operator  and  the  adaptive  crossover  in  Eq  (3.21)  is  to  preserve  individuals  with 
the  best  fitness  and  to  completely  disrupt  individuals  with  the  worst  fitness. 
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3.4.3  Messy  Genetic  Algorithms.  Messy  GAs  (mGAs)  were  originally  developed 
by  Goldberg  [21,  22]  to  overcome  the  problem  of  the  sGA  converging  to  local  optima.  As 
mentioned  in  section  3.3,  for  a  string  length,  m,  there  are  3™  different  possible  schemata. 
Unfortunately,  for  a  population  of  size  n,  only  2™  to  n2™  schemata  are  represented.  This 
problem  is  overcome  by  mGAs.  The  algorithm  is  divided  into  two  phases.  The  first  phase  is 
known  as  a  tournament  and  the  second  phase  is  similar  to  traditional  GAs.  The  tournament 
generates  a  large  number  of  schemata  of  various  sizes  and  then  reduces  this  population  of 
schemata  down  to  a  manageable  size  for  the  second  phase.  The  second  phase  takes  the  reduced 
population  and  uses  genetic  operators  splice,  cut,  and  mutation  to  achieve  an  optimal  solution. 
The  operators  splice  and  cut  replace  the  crossover  operator.  See  [21,  34]  for  a  description 
of  the  splice  and  cut  operators.  Goldberg  states  that  the  mGA  was  able  to  find  the  optimal 
solution  to  a  difficult  problem,  where  the  sGA  only  found  the  optimal  solution  25  per  cent  of 
the  time  [21].  Developments  and  improvements  of  the  mGA  have  been  accomplished  here  at 
AFIT  [33, 16]. 

3.5  Software  Review 

This  section  provides  a  survey  of  available  software  that  was  considered  for  this  thesis. 
For  a  complete  review  of  existing  software  available,  the  reader  is  referred  to  [17].  The  GA 
software  available  can  be  divided  into  three  main  categories:  application  specific,  algorithm 
specific,  and  general  purpose  toolkits.  There  are  no  application  specific  software  available  for 
this  research  topic.  The  next  two  sections  will  discuss  software  that  was  investigated  in  the 
two  remaining  categories. 

3.5.1  Algorithm  Specific.  For  algorithm  specific  software  a  source  code  is  provided 
and  the  user  is  able  to  make  alterations  to  the  code.  Typically  the  code  is  in  a  higher  level 
language  like  ‘C’  and  the  user  interfaces  are  rudimentary.  The  following  is  a  list  of  the  software 
that  was  reviewed: 
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GENESIS  (GENEtic  Search  Implementations  System)  was  written  by  John  Grefenstette  and 
has  been  under  development  since  1981.  The  source  code  is  in  ‘C’  and  allows  for  a  high 
degree  of  modifiability  with  large  amount  of  statistical  information.  It  was  primarily 
developed  for  work  in  a  scientific  environment. 

GAUCSD  (Genetic  Algorithm  University  of  California  San  Diego)  was  written  by  Nicol 
Schraudolph.  It  is  based  on  GENESIS  version  4.5,  but  it  provides  a  higher  level 
of  abstraction  for  defining  the  evaluation  function.  It  allows  direct  use  of  most  ‘C’ 
functions  and  it  has  parallel  capabilities. 

3.5.2  General  Purpose  Toolkits.  These  software  systems  provide  toolkits  of  ac¬ 
cessories  that  can  be  used  interchangeably.  The  toolkits  include  different  graphics  interface 
and  a  library  of  genetic  operators.  The  software  is  designed  to  have  a  user  friendly  inter¬ 
face.  The  library  of  genetic  operators  enables  the  user  to  experiment  with  different  operator 
combinations.  Following  is  a  list  of  the  software  reviewed. 

SPLICER  was  developed  by  Software  Technology  Branch  of  Information  Systems  Direc¬ 
torate  at  NASA/Johnson  Space  center  with  support  from  MITRE  Corporation.  It  has  a 
modular  architecture  that  works  using  Xwindows.  It  was  programmed  in  ‘C’  and  has 
graphics  output  capabilities. 

GAME  (Genetic  Algorithm  Manipulation  Environment)  is  an  algorithm  used  by  a  large 
majority  of  the  European  community.  It  is  designed  for  parallel  capabilities  in  the 
‘C+-I-’  programming  language. 

3.6  Implementation 

3.6.1  Software  Selection.  GENESIS  was  selected  as  the  software  package  to  use 
for  this  thesis  for  a  couple  of  reasons.  First,  it  is  a  readily  available  public  domain  software. 
Second,  work  has  been  accomplished  here  at  the  Air  Force  Institute  of  Technology  (AFIT) 
that  has  modified  the  basic  GENESIS  code  into  a  mGA  [33,  16].  With  these  improvements 
any  work  accomplished  in  this  thesis  could  easily  be  implemented  with  the  modified  codes  to 
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improve  the  computation  time.  The  goal  of  the  implementation  is  to  make  it  simple  to  take  a 
controller  design  and  optimize  it  with  as  little  additional  work  as  possible.  For  this  reason  the 
evaluation  of  the  objective  function  will  be  accomplished  in  the  Matlab®  environment.  This 
is  convenient  since  Matlab®  has  a  large  amount  of  controller  evaluation  routines  already 
programmed.  GENESIS  is  coded  in  ‘C’.  To  link  the  two  together,  Matlab®  is  used  as  a 
computational  engine.  The  GA  sends  a  vector  of  the  design  variables  to  Matlab®  and  a  user 
written  m-file  calculates  the  objective  function  and  passes  the  scalar  value  back  to  GENESIS. 

A  number  of  the  GENESIS  files  were  modified  to  facilitate  this  passing  of  variables 
between  programs.  The  modifications  required  to  run  Matlab®  as  a  computational  engine 
are  documented  in  the  Matlab®  user’s  manual  [31].  Both  the  evaluation  files,  the  ‘C’  code 
and  the  m-file,  are  listed  in  Appendix  B. 
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IV.  Sample  Control  Example 

This  chapter  presents  a  simple  control  problem  for  which  a  gain  schedule  is  designed 
and  optimized  using  GAs.  The  results  of  this  design  clearly  demonstrate  that  the  gain  schedule 
design  process  can  be  optimized.  Further,  for  comparison,  the  gain  schedule  is  optimized  with 
calculus  based  algorithms  and  GAs.  Since,  this  is  a  simple  optimization  problem,  the  calculus 
based  methods  do  show  good  results,  but  they  are  useless  on  the  full  optimization  problem 
shown  in  section  4.4.  The  problem  is  defined  in  section  4.1  and  the  optimization  results  are 
presented  in  sections  4.2, 4.3, 4.4,  and  4.5. 

4 .1  Problem  Statement  and  Analysis 

The  sample  problem  considered  is  a  SISO  third-order  linear  parameter  varying  plant 
with  proportional  feedback  shown  in  Fig.  4.1.  The  scheduling  variable  c  €  (0, 10)  determines 
the  plant  dynamics.  The  range  of  c  is  divided  into  several  intervals  Ii,i  =  1 ...  n,  and  a  gain, 
ki,  is  selected  for  each.  The  controller  gain  is  simply  the  gain  corresponding  to  the  interval 
the  current  value  of  c  is  in.  The  objective  is  to  find  the  set  of  gains  that  minimize  the  average 
deviation  of  the  actual  dominant  closed  loop  poles  from  desired  dominant  pole  locations  over 
the  entire  range  of  c  .  For  efficiency,  the  objective  function  was  approximated  by  evaluating 
the  deviation  at  500  equally  spaced  values  of  c  and  adding  them  together.  The  desired  pole 
locations  are  —2  ±  2j.  Therefore,  the  objective  function  is: 

500 

Ji  {max  [Re(Ai)]  +  2}^  -|-  {max  [Im(Aj)]  —  2}^  (4.1) 

2  =  1 

where  A,  are  the  closed  loop  poles  of  the  system  at  a  given  value  of  c. 

There  are  four  increasingly  difficult  variations  on  this  basic  problem.  First,  the  number  of 
intervals  and  their  sizes  are  chosen  and  the  gains  for  each  interval  are  found.  Results  of  this 
optimization  are  in  section  4.2.  Second,  only  the  number  of  intervals  is  chosen  and  the  interval 
size  and  gains  are  optimized.  These  results  are  in  section  4.3.  Third,  all  three  parameters,  the 
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number  of  intervals,  the  size  of  each  interval,  and  the  gain  for  each  interval ,  are  optimized. 
Results  for  this  are  in  section  4.4.  Finally,  the  gains  are  allowed  to  be  a  linear  function  of  the 
scheduling  variable  and  the  coefficients  of  the  function  are  optimized.  These  results  are  in 
section  4.5. 

In  summary,  the  interval  size  and  number  can  be  specified  where  the  gain  of  each  interval 
is  a  design  variable,  or  just  the  number  of  intervals  can  be  specified  and  both  the  interval 
endpoints  and  gains  are  design  variables,  or  the  number  of  intervals  and  their  respective 
endpoints  and  gains  can  be  design  variables  and  trade  off  complexity  (number  of  intervals) 
for  accuracy  of  the  closed  loop  poles. 


> - fcm 

1 

8^+1  0S^(244C)S+6c 

Rant 


Controller 


Figure  4.1  Example  System  Block  Diagram 


4.2  Fixed  Interval  Results 

Initially,  five  equal  intervals  are  chosen,  (0,2),  [2,4),  [4,6),  [6,8),  and  [8,10),  and  the 
optimal  gains  for  each  interval  are  found.  This  optimization  is  done  using  a  simple  GA  with 
the  following  parameters,  Pc  =  0.95,  =  0.01,  popsize  =  50.  Each  design  variable  had  a 

gene  of  length  14  for  a  precision  of  two  decimal  places.  This  resulted  in  a  total  string  length 
of  100.  For  comparison,  the  same  optimization  is  done  with  a  Broydon-Fletcher-Goldfarb- 
Shanno  (BFGS)  Quasi-Newton  method  with  a  mixed  quadratic  and  cubic  line  search  method. 
There  are  five  design  variables  to  be  optimized,  the  gain  for  each  interval,  li,  i  =  1 ...  5  . 
The  results  of  both  the  GA  and  the  BFGS  method  are  shown  in  Table  4.1.  A  mapping  of  the 
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Table  4. 1  Fixed  Interval  Results 


GA 

BFGS 

#  Iterations 

182 

15 

#  Function  Evals 

8323 

128 

Interval  1  Gain 

31.61 

31.61 

Interval  2  Gain 

22.12 

22.12 

Interval  3  Gain 

13.02 

13.02 

Interval  4  Gain 

4.40 

4.41 

Interval  5  Gain 

-3.62 

-3.59 

Obj.  Function 

91.20 

91.20 

Interval  4  Optimization 


Figure  4.2  Exhaustive  Interval  Optimization 


solution  space  of  one  of  the  intervals  is  shown  in  Fig.  4.2.  For  each  fixed  interval  the  solution 
is  very  similar,  with  one  global  optimum  and  one  local  optimum. 

The  GA  required  a  significantly  greater  number  of  function  evaluations  than  the  BFGS 
method.  However,  this  was  expected  since  the  objective  function  is  smooth  and  the  derivative 
is  well-behaved  for  this  example.  Using  a  GA  at  this  level  is  clearly  overkill,  but  it  is  necessary 
to  verify  that  the  expected  results  are  obtained. 

4.2.1  Parameter  Study.  A  restricted  experimental  study  was  performed  to  find 
the  best  combination  of  GA  control  parameters,  Pc  and  Pm-  Twelve  different  combinations 
of  these  parameters  were  tried.  Each  experiment  had  a  population  size  of  50  and  was  run 
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for  10,000  function  evaluations.  Each  experiment  v^^as  performed  ten  times  with  a  different 
initial  population.  The  results  of  each  experiment  are  listed  in  Table  4.2.  Since  GAs  are  not 
efficient  at  obtaining  a  solution  with  high  precision,  each  experiment  measured  the  number 
of  generations  and  function  evaluations  required  for  the  GA  to  converge  to  within  1  per  cent 
of  the  optimal  solution.  Fig.  4.3  shows  the  average  values  and  their  respective  standard 
deviations  for  each  experiment. 

A  number  of  interesting  trends  are  prevalent  from  these  graphs.  First,  as  the  crossover 
probability  is  reduced,  the  number  of  function  evaluations  per  generation  decreases.  This  is  a 
result  of  the  GA  only  evaluating  the  individuals  of  the  population  that  are  new.  The  second 
interesting  trend  is  that  experiment  number  1  had  the  lowest  number  standard  deviation  for 
both  the  number  of  generations  and  the  number  of  function  evaluations,  7.71  and  365  respec¬ 
tively.  However,  the  next  smallest  standard  deviation  values,  13.01  and  561,  respectively,  are 
significantly  higher.  Additionally,  these  values  are  not  from  the  same  experiment.  Since  the 
parameter  combination  of  experiment  number  1  showed  the  greatest  degree  of  consistency 
in  converging  to  a  solution,  these  parameters  are  used  in  the  remainder  of  the  optimization 
processes. 

4.3  Variable  Interval  Results 

The  same  design  problem  is  now  repeated  with  the  interval  end  points  as  additional 
design  variables.  There  are  now  a  total  of  nine  design  variables.  For  comparison,  the  same 
problem  is  optimized  using  a  sequential  quadratic  programming  (SQP)  method.  For  both  the 
GA  and  the  SQP  method  the  interval  end  points  are  constrained  to  be  within  the  range  of 
(0, 10)  and  the  gains  are  constrained  to  be  within  the  range  of  (—50, 50).  The  results  from 
these  optimizations  are  shown  in  Table  4.3. 

In  Table  4.3,  the  first  column  of  the  SQP  category  shows  the  results  of  the  SQP 
optimization  when  given  the  same  initial  condition  as  the  fixed  interval  optimization.  The  SQP 
method  stalled  at  one  solution  point  for  numerous  iterations  trying  to  find  the  correct  direction 
to  proceed.  In  addition,  when  evaluating  the  derivative  of  the  function,  the  derivative  matrix 
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Table  4.2  Parameter  Study  Results 


Trial  No. 

Pc 

Pm 

Average 
Number  of 
Generation 

Average 
Number  of 
Function 
Evaluations 

Average  No. 

Function 
Evaluations 
per  Generation 

1 

0.95 

0.010 

47.40 

2227 

46.97 

2 

0.85 

0.010 

64.40 

2764 

42.99 

3 

0.75 

0.010 

53.70 

4 

0.65 

0.010 

60.70 

2202 

36.36 

5 

0.95 

0.005 

55.40 

2558 

46.37 

6 

0.85 

0.005 

56.70 

2356 

41.78 

7 

0.75 

0.005 

65.20 

38.86 

8 

0.65 

0.005 

33.68 

9 

0.95 

0.001 

45.76 

10 

0.85 

0.001 

122.10 

4565 

38.48 

11 

0.75 

0.001 

46.00 

1684 

36.96 

12 

0.65 

0.001 

116.40 

3384 

29.94 

Table  4.3  Variable  Interval  Results 


SQP 


Tri£ri  No 


(a)  Generations 


Trial  No 


(c)  Functions  Evaluations  per  Generation 


Figure  4.3  Parameter  Study  Results 
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Rxed  Interval  V.  VariaUe  Interval  Distance  to  Desired  Poles 


Figure  4.4  Variable  Interval  vs.  Fixed  Interval  Pole  Placement  Error 

became  ill  conditioned  and  nearly  singular.  The  resulting  solution  of  this  SQP  optimization 
is  a  local  minimum.  The  SQP  method  was  restarted  with  the  fixed  interval  optimal  solution 
as  its  initial  conditions.  Although  the  same  numerical  difficulties  were  encountered  again,  the 
SQP  algorithms  did  eventually  find  an  optimal  solution.  However,  the  numerical  inaccuracies 
encountered  during  the  optimization  make  the  validity  of  the  solution  questionable.  The  GA 
did  not  have  these  problems  because  it  is  not  dependent  upon  the  local  behavior  of  the  objective 
function  or  its  derivative. 

The  second  derivative  of  the  objective  function  was  evaluated  to  gain  insight  into  the 
shape  of  the  solution  space.  Since  an  analytical  derivative  is  not  practical,  a  second  order 
finite  difference  was  used  to  calculate  the  diagonal  elements  of  the  Hessian.  At  the  optimal 
values  of  the  fixed  interval  solution  found  from  the  BFGS  method,  the  diagonal  elements  are 
[—6  X  10^  —3  X  10^  —1.5  X  10^  —1.5  X  10^  0  0  0  0  0].  The  first  four  terms  are  the 
derivatives  with  respect  to  the  interval  end  points,  and  the  last  five  terms  are  the  derivatives 
with  respect  to  the  gains  of  each  interval.  From  these  results  we  conclude  that  the  solution 
space  is  a  long,  steep  n-dimensional  trough. 
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Figure  4.5  Locus  of  Roots  of  the  Closed  Loop  System 

Fig.  4.4  shows  how  the  pole  placement  error  varies  with  c  for  both  the  variable  interval 
solutions  and  the  fixed  interval  solutions.  Since  the  objective  of  the  optimization  is  to  minimize 
the  distance  from  the  desired  dominant  closed-loop  pole  locations,  the  time  response  of  the 
system  at  various  values  of  the  scheduled  variable  should  be  very  similar.  Fig.  4.5  shows 
a  root  locus  of  the  closed  loop  system  as  c  varies  from  (0 . . .  10).  Fig.  4.6  shows  the  step 
response  for  different  values  of  the  scheduled  variable  from  (0 ...  10)  in  increments  of  0.5. 
As  expected,  the  time  response  in  terms  of  overshoot,  rise  time,  and  settling  time  are  all  quite 
similar.^ 

4.4  Variable  Number  of  Intervals  Results 

For  this  optimization  the  objective  function  is  changed  to  penalize  the  number  of 
intervals.  The  parameters  of  the  sample  problem  are  the  same  as  those  for  the  variable  interval 
case  in  section  4.3,  except  that  the  number  of  intervals  is  allowed  to  vary  between  two  and 
nine.  For  this  example,  the  number  of  design  variables  varies.  If  there  are  only  two  intervals, 
there  are  three  design  variables  (one  interval  end  point,  and  two  controller  gains).  Similarly,  if 
there  are  nine  intervals,  there  are  seventeen  design  variables.  Hence,  if  N  intervals  are  chosen, 
there  are  2N-1  design  variables.  Two  different  objective  functions  are  used  to  observe  the 

^The  steady  state  error  shown  could  be  corrected  with  a  PI  controller. 
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Figure  4.6  Time  Response  for  Various  Values  of  Scheduled  Variable 

effects  of  the  penalty  functions;  Eq  (4.2)  uses  a  quadratic  penalty  function,  and  Eq  (4.3)  uses 
a  linear  penalty  function. 

500 

J2  =  {inax  [Re(Ai)]  +  2}^  +  {max  [Im(Ai)]  -  2}^  +  {N  -  1)^  (4.2) 

2  =  1 

500 

J3  =  ^  {max  [Re(Ai)]  +  2}^  +  {max  [Im(A,)]  -  2}^  +  (iV  -  1)  (4.3) 

2  =  1 

where  A^  is  the  closed  loop  poles  of  the  system  at  a  given  value  of  c. 

Both  equations  were  optimized  using  the  GA  with  the  same  parameters  presented  in 
section  4.2.  The  results  of  these  optimizations  are  shown  in  Table  4.4.  A  calculus  based 
algorithm  was  not  able  to  optimize  this  problem  because  of  the  varying  number  of  design 
variables. 

These  results  show  that  the  initial  arbitrary  use  of  five  intervals  is  too  large  given  these 
objective  functions.  Of  course,  this  result  depends  on  the  relative  weighting  between  the 
two  parts  of  the  objective.  A  more  general  objective  function  would  include  a  multiplicative 
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Table  4.4  Variable  Number  of  Intervals  Optimization  Results 


J2 

#  Iterations 

603 

853 

#  Function  Evals 

18100 

27520 

#  of  Intervals 

4 

4 

Endpoint  1 

2.01 

2.24 

Endpoint  2 

4.31 

4.57 

Endpoint  3 

6.86 

7.11 

Interval  1  Gain 

31.60 

31.33 

Interval  2  Gain 

21.50 

20.64 

Interval  3  Gain 

10.68 

9.86 

Interval  4  Gain 

-0.89 

-1.51 

Obj.  Function 

105.3 

98.6 

weighting  term  in  the  penalty  function.  It  is  interesting  to  note  that  there  was  not  a  significant 
increase  in  computation  time  to  do  this  analysis  and  optimization  of  the  number  of  intervals. 

4.5  Linear  Interpolation 

Since  the  trend  in  the  optimized  gains  is  linear,  a  final  optimization  is  accomplished. 
The  controller  gain  is  chosen  to  be  a  linear  function  of  the  scheduling  variable  c. 

K  =  cti  c  +  ao  (4.4) 

The  coefficients  of  the  linear  function,  ai  and  a2,  are  the  design  variables.  The  BFGS  calculus 
method  used  in  section  4.2  was  also  used  here.  The  GA  parameters  are  Pc  =  0.95,  =  -001, 

and  popsize  =  50.  Table  4.5  compares  the  results  of  the  GA  and  the  BFGS  method.  A  plot  of 
the  root  locus  is  shown  in  Fig.  4.7.  Since,  this  is  a  simple  system,  the  gain  varies  linearly  with 
the  scheduling  variable  c.  Because  of  this  direct  relation  the  calculus  based  method  quickly 
finds  the  optimal  set  of  coefficients.  The  probabilistic  nature  of  the  GA  requires  more  function 
evaluations  to  hone  in  on  the  optimal  solution. 
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4.6  Summary 


This  work  has  shown  some  distinct  advantages  to  using  a  GA  for  gain  schedule  opti¬ 
mization.  For  the  simple  fixed  interval  case,  the  GA  reached  the  same  solution  as  a  calculus 
based  method.  However,  when  the  interval  size  was  allowed  to  vary  the  calculus  based 
method  stalled  and  its  result  was  dependent  on  the  initial  conditions.  More  importantly,  the 
GA  allowed  further  trade-offs  between  the  basic  objective  function  and  the  number  of  intervals 
with  practically  no  increase  in  computational  effort.  Considering  these  key  points,  the  GA 
definitely  shows  promise  for  application  to  real  world  gain  scheduling  problems. 
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Table  4.5  Linear  Scheduling  Results 


BFGS 

#  Iterations 

48 

7 

#  Function  Evals 

1954 

38 

«i 

-4.377 

ao 

34.946 

35.1139 

Obj.  Function 

82.31 

82.30 

Figure  4.7  Locus  of  Roots  of  the  Closed  Loop  System 
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V.  Flight  Envelope  Design  Example 

The  previous  chapter  demonstrated  the  applicability  of  gain  scheduling  optimization  on 
a  simple  SISO  system.  Now,  an  optimization  of  a  real  world  system  is  completed.  A  full 
envelope  controller  for  an  F-18  Supermanueverable  fighter  in  Fig.  5.1  was  designed  in  [4] 
using  a  simple  gain  scheduling  technique  of  linear  interpolation.  This  chapter  optimizes  this 
gain  schedule  using  the  methods  in  Chapter  IV.  First,  the  linearized  equations  of  motion  are 
derived  in  section  5.1.  Next,  the  form  of  the  controller  and  the  gain  schedule  developed  in  [4] 
are  presented  in  section  5.2  as  a  baseline  for  optimization.  Afterwards,  a  measure  of  relative 
error  is  presented  and  the  optimization  function  is  developed  in  section  5.3. 

5.1  Equations  of  Motion 

The  nonlinear  equations  of  motion  of  an  aircraft  are  shown  in  Eq  (5.1)  and  Eq  (5.2).  Eq 
(5.1)  describes  the  forces  along  the  aircraft  body  axes,  and  Eq  (5.2)  describes  the  rotational 
forces  about  the  same  axes. 

=  m{U  ->r  WQ  -  V R  +  gsme) 

Fy  =  m{V  +  UR- WP -gcosQsm^)  (5.1) 

F,  =  m{W  +  VP-UQ-  cos  0COS  $) 

L  =  PIx  RIxz  +  QR{Iz  ~  ly)  ~  PQIxz 
M  =  QIy  +  PR{Ix  -  Q  -  RHx,  +  P^Ixz  (5.2) 

N  =  RR  +  P Ixz  T  PQ{Iy  —  Iz)  +  QRIxZ 

Fx,  Fy,  and  F^  are  the  external  forces  and  L,  M,  and  N,  are  the  external  moments. 

U,  V,  W  are  the  translational  velocities  along  the  X,  Y,  and  Z  body  axes  directions 
respectively.  P,  Q,  and  R  are  the  rotational  velocities  about  the  X,  Y,  and  Z  axes, 
respectively.  E,  I y,  I z,  and  Ez  are  the  moments  of  inertia,  g  is  gravity,  and  m  is 
the  aircraft  mass. 
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Figure  5.1  F-18  Aircraft 


Since  these  forces  and  moments  are  in  the  body  axis  system,  Eq  (5.3)  relates  the  aircraft 
orientation  to  the  Earth  through  the  Euler  angles  0,  and 


Q  =  Q  cos  $  —  J?  sin  # 


$  =  P  4-  (5  tan  0  sin  $  +  i?  tan  0  cos  $ 

•  B  cos  $  Q  sin  $ 
qr  = - - 

cos  0  cos  0 


(5.3) 


where  0  is  the  pitch  angle,  $  is  the  roll  angle,  and  ’S'  is  the  yaw  angle. 

These  nonlinear  equations  are  linearized  assuming  the  following: 

•  The  aircraft  is  a  rigid  body, 

•  The  aircraft  mass  is  constant, 

•  The  aircraft  is  trimmed  to  an  equilibrium  condition  such  that  all  accelera¬ 
tions  are  zero, 

•  The  aircraft  is  in  straight  and  level  flight, 

•  The  aircraft  dynamics  can  be  decoupled  into  longitudinal  and  lateral/directional 
components, 

•  The  linear  and  rotational  velocities,  and  the  external  forces  and  moments 
can  be  considered  perturbed  from  equilibrium  values. 
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•  The  product  of  small  perturbed  values  are  small  and  negligible,  and 

•  The  angles  between  the  the  equilibrium  and  the  disturbed  value  are  consid¬ 
ered  small. 

Additionally,  since  the  short  period  mode  dominates  the  aircraft  response  to  pilot  inputs, 
the  longitudinal  equations  are  simplified  as  shown  in  Eq  (5.4).  Further  explanation  of  these 
assumptions  and  simplifications  can  be  found  in  [4,  8]. 


a 

Za 

' 

a 

+ 

'  Zs^ 

^SpTV 

Se 

q 

M, 

q 

Ms^ 

■^Spxv 

SpTV 

where  o.  is  the  angle  of  attack,  q  is  the  linearized  pitch  rate,  is  the  elevator 
deflection,  is  the  pitch  thrust  vectoring,  Z  and  M  are  longitudinal  stability 
derivatives. 

Through  the  use  of  a  nonlinear  control  selector,  developed  in  [4],  the  thrust  vectoring 
commands  are  only  used  when  the  aerodynamic  forces  are  not  sufficient  to  complete  the  com¬ 
manded  maneuver.  Furthermore,  the  use  of  the  control  selector  enables  the  conunanded  inputs 
to  be  transformed  into  generalized  command  inputs:  pitch  acceleration,  pc,  yaw  acceleration, 
qc,  and  roll  acceleration,  rV 

From  this  mapping,  the  longitudinal  linear  design  model  is  that  shown  in  Eq  (5.5). 
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where  Biong  = 

1 

Ciong  — 

0 

1 

(5.5) 

(5.6) 


where  Along  is  the  longitudinal  plant  matrix,  Biong  is  the  longitudinal  input  matrix, 
and  Ciong  is  the  longitudinal  output  matrix. 
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Figure  5.2  Outer  Control  Loop 
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Figure  5.3  Inner  Equalization  Loop 


5.2  Controller  Design 

The  control  of  the  F-18  longitudinal  dynamics  is  decomposed  into  an  inner  and  outer 
loop  (See  Fig.  5.3  and  Fig.  5.2).  The  objective  of  the  inner  loop  is  to  equalize  the  closed 
loop  frequency  response  over  the  flight  envelope  so  that  the  outer  loop  controller  does  not 
need  to  be  scheduled.  Therefore,  this  investigation  focuses  on  the  inner  loop  controller.  The 
original  design  selected  twelve  flight  conditions,  shown  in  Fig.  5.4,  at  which  to  design  an 
inner  loop  controller.  At  each  design  point  a  minimal  i^oo  controller  was  designed  [52].  Then, 
using  engineering  judgement,  a  central  controller  was  selected  as  Mach  0.95  and  an  altitude 
of  20000  feet.  The  closed  inner  loop  of  this  central  controller,  Po,  is  used  to  design  the  outer 
loop  controller,  Kout- 


5-4 


Altitude 

(ft) 


xlO’ 


Figure  5.4  F-18  Flight  Envelope 


The  inner  loop  controller,  shown  in  Fig.  5.5,  has  three  parameters  that  are  scheduled 
with  dynamic  pressure,  q. 

F-GN  Kf-GM 

K  - _ - _ 

-N  -M 

where  F  -  [-40]  Kf  =  [ii]  G=  [.0247] 

N  =  N{q)  M  =  [  Mr{q)  M^iq)  ]  (5-7) 

The  parameters  N  and  M  are  scheduled  as  a  linear  function  of  dynamic  pressure.  The  schedules 
for  these  parameters  are  shown  in  Eq  (5.8)  and  graphically  in  Fig.  5.6. 

N{q)  =  -0.312^  +  461 

Ml  (g)  = -0.0589  +  50.5  (5.8) 

M2(9)  =  -O.OO69  +  8.II 


Keq 


Figure  5.5  Inner  Loop  Controller 
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The  schedules  were  developed  by  plotting  the  H^o  scheduled  controller  parameters  as 
a  function  of  dynamic  pressure.  A  linear  fit  of  the  parameters  was  performed  using  a  least 
square  error  method.  These  controller  schedules  are  the  baseline  design  for  comparison  with 
the  ‘optimized’  schedules.  The  goal  of  the  optimization  is  to  improve  this  schedule  in  terms 
of  relative  error,  which  is  described  in  the  next  section. 


1000 


Figure  5.6  Baseline  Schedule  of  Design  Parameters 


5.3  Relative  Error 

Relative  error,  defined  by  Eq  (5.9),  is  used  to  measure  the  equalization  of  the  closed 
inner-loop. 

Am  =  {P-  Po)Po^  (5.9) 


Po  is  the  closed  inner  loop  used  for  outer  loop  controller  design,  and  P  is  a  closed 
inner  loop  at  another  flight  condition. 

Safonov  and  Chiang’s  Robustness  Theorem  [42]  provides  a  weak  sufficient  condition  for 
stability  using  relative  error.  The  authors  in  [4]  succinctly  stated  this  sufficient  condition. 

If  a{Km)  <  liox  uj  <  Ur,  then  the  closed  loop  system  will  be  stable  provided 
that  the  control  bandwidth,  Ub,  is  less  than  Ur. 
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Therefore,  any  outer  loop  controller  designed  for  Pq  will  be  stable  for  any  inner  loop  plant, 
P,  provided  the  relative  error,  A^,  is  sufficiently  small. 

A  summary  of  the  proof  was  presented  in  [4]  and  is  reproduced  here  for  completeness. 

Consider  the  system  with  equalized  plant,  P,  outer  loop  controller,  Koi,  and  rela¬ 
tive  error,  A^.  The  bandwidth  of  the  control  system  is  defined  as  the  frequency 
range  where  the  loop  transfer  function  is  big.  That  is: 

a{PKoi)  >1  yu}<LOb  (5.10) 

A  sufficient  condition  for  stability  is 

a{AmMPKoi{I  +  PKoi)-^)  <  1  (5.11) 

It  follows  from  Eq  (5.1 1)  that  for  to  <  Ub, 

a{PKM  +  PKoi)'")-^  (5.12) 

So,  for  frequencies  where  the  loop  transfer  function  gain  is  big,  a  sufficient 
condition  for  stability  is 

H^m)  <  1  (5.13) 

Relative  error  is  a  frequency  domain  measure  of  the  differences  between  closed  loop 
system  responses  at  various  flight  conditions.  The  main  goal  of  the  inner-loop  of  this  design 
is  to  achieve  a  relative  error  less  than  one  for  all  flight  conditions.  The  logical  optimization 
objective  is  to  minimize  this  error  for  all  flight  conditions  and  thereby  achieve  a  more  equalized 
inner-loop.  Therefore,  the  objective  function  for  this  optimization  is 

min  J  =  ^d-(Am,.)  (5.14) 

i=l 

where  is  defined  in  Eq  (5.9),  and  n  is  the  number  of  flight  conditions  to  be 
evaluated.  ^ 

The  flight  conditions  used  in  evaluating  the  objective  function  are  the  original  twelve  used  for 
design  in  [4]  plus  an  additional  eight  flight  conditions  chosen  to  more  uniformly  represent  the 

^For  n  flight  conditions  there  are  only  {n  —  1)  nonzero  relative  errors. 
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flight  envelope  as  a  function  of  dynamic  pressure.  The  original  twelve  flight  conditions  are 
shown  in  Fig.  5.4.  All  of  these  flight  conditions  and  their  plant  matrices  are  listed  in  Appendix 
A.  The  goal  of  this  optimization  is  that  by  reducing  the  relative  error  of  the  inner  loop  the 
time  response  of  the  system  will  become  more  uniform  throughout  the  operating  envelope. 

5.4  Optimization  Results 

This  section  presents  the  optimization  results  of  the  F-18  longitudinal  inner  equaliza¬ 
tion  loop  controller.  Multiple  optimizations  were  performed  with  various  alterations  on  the 
basic  objective  function  given  in  Eq  (5.14).  The  first  section  provides  on  overview  of  the 
optimization  cases  considered  and  the  respective  sections  in  which  the  results  are  presented. 
Next,  section  5.4.11  compares  and  summarizes  the  optimization  results,  and  section  5.4.12 
demonstrates  the  responses  of  off  design  flight  conditions  using  the  optimized  gain  schedule. 

5.4.1  Overview.  The  optimization  process  is  performed  with  a  family  of  twenty 
flight  conditions  representing  the  flight  envelope  shown  in  Fig.  5.7.  The  linearized  plants  are 
listed  in  Appendix  A.  The  optimization  results  are  based  on  four  objective  functions 

20 

min  Ji  =  (5.15) 

i=l 

Note:  Since  the  relative  error  is  a  measure  of  relative  difference,  one  of  the 

relative  error  values  is  zero. 

20 

min  J2  =  XI  +  f{N)  (5.16) 

where  f{N)  is  a  functional  weighting  on  the  number  of  intervals. 

20 

min  J3  =  +  ^«(^) 

i=l 

where  ei{t)  is  the  error  between  the  time  response  of  the  chosen  central  flight 

condition  and  the  time  response  of  the  flight  condition 
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20  ^  1 

max  J4  =  +  -ei{t)  (5.18) 

i=i  ^ 

Eight  optimization  cases  are  considered  and  are  summarized  in  Table  5.1.  A  brief 
explanation  of  each  optimization  is  presented  in  the  following  list: 

Case  I  The  scheduled  variables  are  piecewise  constant  functions  of  the  scheduling  variable. 
Four  intervals  of  the  scheduling  variable  were  arbitrary  chosen  to  be  of  equal  size: 
(0, 250],  (250, 500],  (500, 750],  (750, 1000).  Ji  is  the  objective  function.  The  central 
flight  condition,  Pq,  was  chosen  as  Mach  0.95  and  an  altitude  of  20,000  feet  [4]. 

Case  II  This  case  is  the  same  as  Case  I  except  that  the  size  of  the  interval  is  allowed  to  vary 
to  reduce  the  objective  function  further. 

Case  III  This  case  is  the  same  as  Case  II  except  that  the  number  of  intervals  is  also  a  design 
variable.  J2  is  the  objective  function. 
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Case 

Objective  Function 

Results  Section 

0 

Baseline 

5.4.2 

I 

Ji 

5.4.3 

II 

Ji 

5.4.4 

m 

J2 

5.4.5 

IV 

J2 

5.4.6 

V 

J3 

5.4.7 

VI 

Ji 

5.4.8 

vn 

J3 

5.4.9 

vm 

J4 

5.4.10 

Table  5.1  Optimization  Summary 


Case  IV  This  case  is  the  same  as  Case  III  except  the  central  flight  condition  is  allowed  to 
vary  within  the  family  of  plants  to  further  reduce  the  objective  function,  J2. 

Case  V  This  case  uses  objective  function  J3.  The  central  flight  condition,  the  number  of 
intervals,  and  the  size  of  each  interval  are  allowed  to  vary. 

Case  VI  For  this  case  the  scheduling  variables  are  piecewise  linear  functions  of  dynamic 
pressure.  The  objective  function  of  the  optimization  is  J\.  The  central  flight  condition, 
number  of  intervals,  and  the  size  of  the  intervals  are  allowed  to  vary. 

Case  VII  This  case  is  the  same  as  Case  VI  except  that  the  objective  function  is  J3. 

Case  VIII  This  case  is  the  same  as  Case  VI  except  that  the  objective  function  is  J4. 

5.4.2  Baseline  Design  Results.  The  baseline  design  developed  in  [4],  was  not 
optimized.  Relative  error  was  only  used  as  a  weak  sufficient  condition  for  stability  of  the 
controller  throughout  the  flight  envelope.  The  scheduling  variables  were  chosen  as  described 
in  section  5.2.  The  open  loop  frequency  responses  of  the  inner  loop  are  shown  in  Fig.  5.8  and 
Fig.  5.9.  The  flight  condition  with  the  lowest  dynamic  pressure  has  the  highest  low  frequency 
gain  of  aU  the  flight  conditions.  The  general  trend  is  that  as  the  dynamic  pressure  increases 
the  low  frequency  gain  decreases. 
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Figure  5.8  Case  0;  qc  to  a  Open-Loop  Dynamics 


Figure  5.9  Case  0:  qc  to  q  Open-Loop  Dynamics 
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Figure  5.10  Case  0:  qc  to  a  Closed-Loop  Dynamics 


The  closed  loop  dynamics  of  the  inner  loop  are  shown  in  Fig.  5.10  and  Fig.  5.11.  The 
two  most  important  graphs  for  demonstrating  the  optimization  results  are  Fig.  5.12  and  Fig. 
5.13.  These  graphs  depict  the  relative  error  of  the  closed  inner  equalization  loop  and  the  time 
response  to  a  step  pitch  acceleration  command  of  the  closed  outer  loop,  respectively.  Note  that 
the  time  response  has  two  distinct  groupings.  This  is  because  there  are  actually  two  outer-loop 
controllers,  one  for  low  dynamic  pressure  {q  <  200  psf)  and  one  for  high  dynamic  pressure 
(q  <  200  psf).  This  is  a  consequence  of  two  desired  pitch  responses  in  different  regions  of  the 
flight  envelope.  At  higher  velocities,  the  pilot  likes  to  feel  a  faster  time  response.  Also  note 
that  the  relative  error  for  most  flight  conditions  has  a  maximum  singular  value  around  0.4  to 
0.5. 

The  baseline  schedules  of  the  parameters  N  and  M  are  shown  in  Fig.  5.14.  The  five 
graphs.  Fig.  5.10,  Fig.  5.11,  Fig.  5.12,  Fig.  5.13,  and  Fig.  5.14  are  the  baseline  for 
comparison  with  the  optimization  results. 

The  values  of  the  objective  functions  are  Ji  =  6.37,  J3  =  23.24,  and  J4  =  9.74. 
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Singular  Value  Singular  Value 


Figure  5.1 1  Case  0:  4c  to  q  Closed-Loop  Dynamics 


w  (rad/sec) 


Figure  5.12  Case  0:  Relative  Error 
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Figure  5.13  Case  0:  Time  Response  for  a  Step  Input 


Figure  5.14  Case  0:  Controller  Parameter  Schedule 
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5.4.3  Case  I.  The  design  variables  of  the  GA  optimization  are  the  scheduled 
parameters  defined  in  Eq  (5.8).  The  schedules  are  evaluated  at  gmin  =  0  and  ^max  =  1000 
to  obtain  the  minimum  and  maximum  values  of  the  parameters  N  and  M.  For  a  first  look  at 
optimizing  these  parameters,  the  range  of  the  scheduling  variable,  q,  is  arbitrarily  divided  into 
four  intervals,  (0,250],  (250,500],  (500,750],  and  (750,1000).  The  scheduled  parameters  are 
chosen  to  be  piecewise  constant  values  of  dynamic  pressure  in  each  interval.  The  parameters 
N  and  M  are  optimized  for  each  interval  to  minimize  objective  function  Ji.  There  are  twelve 
design  variables  to  optimize. 

The  optimized  objective  function  value  is  3.09,  half  of  the  baseline  relative  error.  The 
closed  inner-loop  responses  are  shown  in  Fig.  5.15  and  Fig.  5.16.  The  relative  error  and  time 
responses  are  shown  in  Fig.  5.17  and  Fig.  5.18.  The  resulting  scheduled  variables  are  shown 
in  Fig.  5.19. 

As  a  result  of  the  reduction  in  relative  error,  the  overshoot  in  the  time  response  has 
decreased  for  most  of  the  flight  conditions.  However,  for  three  flight  conditions  the  overshoot 
increased  over  the  baseline.  This  can  also  be  seen  in  the  relative  error  graph.  Fig.  5.17, 
where  three  relative  errors  are  increased  over  the  baseline  and  the  remaining  relative  errors  are 
reduced.  Additionally,  the  closed  inner-loop  pitch  response  is  more  uniform  than  the  baseline 
design. 

5.4.4  Case  II.  For  this  optimization  the  interval  size  is  allowed  to  vary.  The  number 
of  intervals  is  still  constrained  to  be  four.  Now  the  number  of  design  variables  is  fifteen,  three 
interval  endpoints  and  twelve  scheduled  variable  values.  Each  interval  endpoint  is  allowed  to 
vary  between  0  and  1000  psf. 

The  objective  function  value  is  2.51  which  is  less  than  Case  I.  The  closed  inner-loop 
responses  are  shown  in  Fig.  5.20  and  Fig.  5.21.  The  relative  error  and  time  responses  are 
shown  in  Fig.  5.22  and  Fig.  5.23.  The  resulting  scheduled  variables  are  shown  in  Fig.  5.24. 
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Figure  5.17  Case  I:  Relative  Error 


Time  (sec) 

Figure  5.18  Case  I:  Time  Response  for  a  Step  Input 
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Figure  5.19  Case  I:  Controller  Parameter  Schedule 

As  in  Case  I,  the  time  response  of  two  flight  conditions  have  overshoots  larger  than 
the  baseline.  The  significance  of  this  optimization  is  that  the  relative  error  is  decreased  by 
allowing  the  interval  size  to  vary. 

5.4.5  Case  III.  For  this  optimization  the  number  of  intervals  is  allowed  to  vary 
between  (2..9).  Additionally,  the  interval  endpoints  are  allowed  to  vary  between  (0..1000). 
Therefore  the  number  of  design  variables  varies  between  (8. .36).  Four  different  penalty 
functions  on  the  number  of  intervals  are  optimized. 

a)/i(7V)  =  iV 
h)h{N)  =  N^ 

c) f3{N)  =  Nl9 

d) MN)  =  {N/9Y 

The  results  of  each  of  these  optimizations  are  in  the  following  subsections. 

5. 4.5.1  Case  Ilia.  Two  intervals  are  chosen  as  with  an  objective  function 
value  of  5.24.  The  closed  inner-loop  responses  are  shown  in  Fig.  5.25  and  Fig.  5.26.  The 
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Figure  5.22  Case  H:  Relative  Error 
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Figure  5.24  Case  11:  Controller  Parameter  Schedule 

relative  error  and  time  responses  are  shown  in  Fig.  5.27  and  Fig.  5.28  and  the  resulting 
scheduled  variables  are  shown  in  Fig.  5.29. 

The  time  responses  show  a  significant  improvement  in  overshoot  when  compared  with 
Cases  I  and  H.  However,  there  are  still  three  time  responses  with  overshoots  larger  than  the 
baseline.  The  rise  time  and  settling  time  are  consistent  with  variations  in  the  flight  condition. 
The  relative  error  is  about  the  same  value  as  Case  I  when  the  penalty  function  value  is 
subtracted  off.  The  significant  finding  is  that  the  same  minimization  of  relative  error  can  be 
achieved  with  half  the  intervals.  Therefore,  the  same  improvements  in  the  equalization  of  the 
inner-loop  can  be  achieved  with  half  the  complexity. 


5.4.52  Case  mb.  Two  intervals  are  again  chosen  as  optimal  with  an  objective 
function  value  of  7.23.  The  closed  inner-loop  responses  are  shown  in  Fig.  5.30  and  Fig.  5.31. 
The  relative  error  and  time  responses  are  shown  in  Fig.  5.32  and  Fig.  5.33  and  the  resulting 
scheduled  variables  are  shown  in  Fig.  5.34. 
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Figure  5.27  Case  Ma:  Relative  Error 


Time  (sec) 

Figure  5.28  Case  Ilia:  Time  Response  for  a  Step  Input 
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Figure  5.29  Case  Ola:  Controller  Parameter  Schedule 

There  is  little  to  no  difference  in  the  closed  loop  responses,  the  relative  error,  and  the 
time  responses  between  this  case  and  Case  nia.  The  only  difference  is  a  small  change  in  the 
scheduled  parameter  N.  The  parameters  Ml  and  M2  are  even  the  same  for  these  two  cases. 
This  seems  to  indicate  that  the  parameter  N  is  not  very  significant  in  the  scheduling  process. 

5.4.53  Case  IIIc.  Two  intervals  are  again  chosen  with  an  objective  function 
value  of  3.45.  The  closed  inner-loop  responses  are  shown  in  Fig.  5.35  and  Fig.  5.36.  The 
relative  error  and  time  responses  are  shown  in  Fig.  5.37  and  Fig.  5.38  and  the  resulting 
scheduled  variables  are  shown  in  Fig.  5.39. 

For  this  case  the  closed  loop  responses,  the  relative  error,  the  time  responses,  and  the 
scheduling  variables  N,  Ml  and  M2  all  show  the  same  trends  as  Case  Dla.  Again,  N  is  slightly 
different. 


5. 4. 5. 4  Case  Illd.  Here,  three  intervals  are  chosen  with  and  objective  function 

value  of  3.24.  Since  the  penalty  for  the  number  of  intervals  is  the  least  for  this  optimization,  a 
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Figure  5.32  Case  nib:  Relative  Error 
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Figure  5.33  Case  Illb:  Time  Response  for  a  Step  Input 
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Figure  5.37  Case  nic:  Relative  Error 
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Figure  5.38  Case  IIIc:  Time  Response  for  a  Step  Input 


Figure  5.39  Case  IIIc:  Controller  Parameter  Schedule 
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Figure  5.40  Case  Hid:  qc  to  a  Closed-Loop  Dynamics 


larger  number  of  intervals  is  selected.  The  closed  inner-loop  responses  are  shown  in  Fig.  5.40 
and  Fig.  5.41.  The  relative  error  and  time  responses  are  shown  in  Fig.  5.42  and  Fig.  5.43  and 
the  resulting  scheduled  variables  are  shown  in  Fig.  5.44. 

The  only  significant  difference  between  this  case  and  Cases  rHa,  b,  and  c  is  the  time 
response.  The  time  response  resulting  from  this  optimization  is  comparable  to  the  time 
response  of  Case  I.  The  time  responses  for  this  case  are  slightly  worse  than  Cases  Ilia,  b,  and 
c. 


5.4.6  Case  IV.  For  this  optimization,  the  central  flight  condition  which  was  chosen 
by  engineering  judgement  for  the  baseline  design  is  allowed  to  vary.  The  number  of  design 
variables  is  increased  by  one;  the  number  of  design  variable  now  varies  between  9  and  37. 
Additionally,  two  different  penalty  functions  on  the  number  of  intervals  are  optimized. 

a)  /i(W)  =  0 

b)  /2(A')  =  Af/100 

The  results  of  each  optimization  are  in  the  following  subsections. 
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Figure  5.42  Case  Hid;  Relative  Error 
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Figure  5.43  Case  IHd:  Time  Response  for  a  Step  Input 
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Figure  5.44  Case  Hid;  Controller  Parameter  Schedule 
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Figure  5.45  Case  IVa:  4c  to  a  Closed-Loop  Dynamics 


5 A. 6.1  Case  IVa.  Five  intervals  are  chosen  with  an  objective  function  value 
of  2.07.  This  can  be  reduced  to  four  intervals  to  eliminate  the  ‘spike’  in  the  schedule.  The 
closed  inner-loop  responses  are  shown  in  Fig.  5.45  and  Fig.  5.46.  The  relative  error  and 
time  responses  are  shown  in  Fig.  5.47  and  Fig.  5.48  and  the  resulting  scheduled  variables  are 
shown  in  Fig.  5.49. 

The  most  significant  result  of  this  optimization  is  that  the  relative  error  is  reduced  by 
allowing  the  central  flight  condition  to  vary.  Note  that  the  worst  relative  error  has  a  value  of 
0.32  and  the  best  relative  error  has  a  value  of  0.0001.  Additionally,  even  though  there  was 
no  penalty  on  the  number  of  intervals,  the  maximum  allowable  number  of  intervals  was  not 
optimal.  The  closed  loop  responses  and  the  time  responses  are  slightly  better  than  those  in 
Case  m. 


5. 4. 6.2  Case  IVb.  Here,  seven  intervals  are  chosen  with  an  objective  function 
value  of  2.26.  However,  this  can  be  reduced  to  six  by  eliminating  the  ‘spike’  in  the  schedules. 
The  closed  inner-loop  responses  are  shown  in  Fig.  5.50  and  Fig.  5.51.  The  relative  error  and 
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Figure  5.47  Case  IVa:  Relative  Error 
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Figure  5.48  Case  IVa:  Time  Response  for  a  Step  Input 
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Figure  5.49  Case  IVa:  Controller  Parameter  Schedule 


Figure  5.50  Case  IVb:  qc  to  a  Closed-Loop  Dynamics 


time  responses  are  shown  in  Fig.  5.52  and  Fig.  5.53  and  the  resulting  scheduled  variables  are 
shown  in  Fig.  5.54.  The  objective  function  value  is  2.26  with  seven  intervals. 

The  results  of  this  optimization  are  similar  to  those  in  Case  IVa.  The  significant 
difference  is  in  the  time  response.  For  the  flight  conditions  evaluated,  the  time  response  for 
all  flight  conditions  was  very  nearly  uniform  even  though  the  resulting  objective  function  is 
not  the  smallest  found  so  far.  This  seems  to  indicate  that  there  is  not  as  close  a  correlation 
between  the  closed  loop  frequency  response  and  the  closed  loop  time  response  as  previously 
thought. 

5.4.7  Case  V.  Since  there  is  an  indication  that  the  closed  loop  frequency  response 
error  and  the  closed  loop  time  response  are  not  directly  related,  this  optimization  directly 
optimizes  both  objectives.  The  time  domain  error  is  measured  as  the  difference  between  the 
central  controller’s  time  response  and  the  time  response  at  all  the  other  flight  conditions.  Time 
time  response  if  evaluated  for  5.0  seconds  at  increments  of  0.1  seconds  to  determine  the  time 
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Figure  5.52  Case  IVb:  Relative  Error 
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Figure  5.54  Case  IVb:  Controller  Parameter  Schedule 
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error  (See  Eq  (5.19)). 


50 


t  =  0,0.1, 0.2,  •■•,5 


(5.19) 


where  yo{t)  is  the  time  response  of  the  closed  outer-loop  at  the  central  flight 
condition  and  yj{t)  is  the  time  response  of  the  closed  outer  loop  at  a  different 
flight  condition. 

The  number  of  design  variables  is  the  same  as  in  Case  IV. 

Five  intervals  are  chosen  with  an  objective  function  value  of  14.22.  The  closed  inner- 
loop  responses  are  shown  in  Fig.  5.55  and  Fig.  5.56.  The  relative  error  and  time  responses 
are  shown  in  Fig.  5.57  and  Fig.  5.58  and  the  resulting  scheduled  variables  are  shown  in  Fig. 
5.59. 

The  closed  loop  frequency  responses  are  similar  to  those  of  the  previous  cases.  The 
relative  error  is  reduced  from  the  baseline  but  is  not  near  the  minimum  relative  error  value 
found  so  far.  However,  the  time  response  is  very  closely  uniform  for  all  of  the  flight  conditions. 

5.4.8  Case  VI.  This  optimization  duplicates  Case  IV  except  that  the  scheduled 
variables  are  now  piecewise  linear  functions  of  the  scheduling  variable.  The  number  of  design 
variables  is  increased  by  three,  since  an  addition  set  of  scheduled  variables  is  needed  for  a 
linear  fit  The  number  of  design  variables  now  varies  between  12  and  40. 

Three  intervals  are  chosen  with  an  objective  function  value  of  1.98,  the  lowest  of  all 
cases.  The  closed  inner-loop  responses  are  shown  in  Fig.  5.60  and  Fig.  5.61.  The  relative 
error  and  time  responses  are  shown  in  Fig.  5.62  and  Fig.  5.63  and  the  resulting  scheduled 
variables  are  shown  in  Fig.  5.64. 

The  significant  difference  between  this  case  and  the  previous  cases  is  that  the  minimum 
relative  error  was  found  sooner  with  fewer  function  evaluations.  This  case  also  found  the 
smallest  relative  error  of  all  the  optimization  cases.  The  result  of  this  minimal  relative  error  is 
evident  in  Fig.  5.60  and  Fig.  5.61  where  there  is  a  tight  grouping  of  the  closed  loop  response. 
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Figure  5.57  Case  V:  Relative  Error 
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Figure  5.58  Case  V:  Time  Response  for  a  Step  Input 
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Figure  5.59  Case  V:  Controller  Parameter  Schedule 

Since  the  schedule  is  better,  piecewise  linear  versus  piecewise  constant,  the  optimization 
results  are  better  than  the  previous  cases. 

5.4.9  Case  VII.  This  optimization  duplicates  Case  V  for  piecewise  linear  scheduled 
variables.  The  number  of  design  variables  now  varies  between  12  and  40. 

Two  intervals  are  chosen  as  optimal  with  no  penalty  function  on  the  number  of  intervals. 
The  objective  function  value  is  12.25,  which  is  less  than  Case  V  results.  The  closed  inner-loop 
responses  are  shown  in  Fig.  5.65  and  Fig.  5.66.  The  relative  error  and  time  responses  are 
shown  in  Fig.  5.67  and  Fig.  5.68  and  the  resulting  scheduled  variables  are  shown  in  Fig.  5.69. 

The  closed  loop  frequency  responses  and  the  relative  error  results  are  similar  to  those 
in  Case  V,  but  the  time  response  is  not  quite  as  uniform  as  the  time  response  in  Case  V.  This 
might  seem  odd  since  the  value  of  the  objective  function  is  less  than  that  of  Case  V.  This  is 
because  each  of  these  cases  is  using  a  different  central  controller. 
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Figure  5.62  Case  VI:  Relative  Error 
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Figure  5.64  Case  VI:  Controller  Parameter  Schedule 


Figure  5.65  Case  VII:  4c  to  a  Closed-Loop  Dynamics 
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Figure  5.66  Case  VII:  qc  to  q  Closed-Loop  Dynamics 


Figure  5.67  Case  VII;  Relative  Error 
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Figure  5.68  Case  VII:  Time  Response  for  a  Step  Input 
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Figure  5.69  Case  VII:  Controller  Parameter  Schedule 
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Figure  5.70  Case  Vni:  qc  to  a  Closed-Loop  Dynamics 


5.4.10  Case  VIII.  For  this  optimization,  objective  function  J4  is  used  to  more 
equally  weight  the  time  domain  error  and  the  frequency  domain  error.  The  design  variables 
are  the  same  as  in  Case  Vn.  Since  the  relative  error  is  less  in  magnitude  than  the  time  domain 
error,  the  time  domain  error  is  scaled  smaller.  Whereas,  objective  function  J3  placed  an 
emphasis  on  the  time  domain  error. 

Two  intervals  are  again  chosen  as  optimal  with  no  penalty  on  the  number  of  intervals. 
The  objective  function  value  is  5.16.  The  closed  inner-loop  responses  are  shown  in  Fig.  5.70 
and  Fig.  5.71.  The  relative  error  and  time  responses  are  shown  in  Fig.  5.72  and  Fig.  5.73  and 
the  resulting  scheduled  variables  are  shown  in  Fig.  5.74. 

The  time  response  is  almost  uniform  over  the  operating  envelope  for  the  points  evaluated. 
Additionally,  the  closed  loop  responses  are  very  uniform  for  frequencies  greater  than  5  radians 
per  second.  Furthermore,  the  control  parameter  schedules  are  very  simple  with  only  two 
intervals  and  piecewise  linear  function  of  dynamic  pressure.  Therefore,  this  optimization  is 
selected  as  the  best.  This  schedule  is  validated  in  section  5.4.12. 
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Figure  5.72  Case  Vni:  Relative  Error 
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Figure  5.73  Case  VIII:  Time  Response  for  a  Step  Input 
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Figure  5.74  Case  Vni:  Controller  Parameter  Schedule 
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5.4.11  Summary  of  Results.  A  summary  of  the  optimization  results  are  shown 
in  Table  5.2.  For  each  case  the  values  of  all  four  objective  functions  are  calculated  where 
possible  for  comparison.  The  best  result  for  each  objective  function  is  emphasized.  No 
single  optimization  is  able  to  optimize  more  than  one  objective  simultaneously.  Each  case 
demonstrates  that  the  gain  schedule  can  be  optimized  with  respect  to  the  specified  objective 
function.  By  changing  the  objective  function  to  meet  the  desired  goal,  an  optimal  gain  schedule 
is  achievable. 

Every  single  optimization  case  tried  is  able  to  improve  on  the  baseline  design.  The 
resulting  optimization  has  a  more  consistent  time  response  for  the  flight  conditions  evaluated. 
Additionally,  another  benefit  of  the  optimization  is  that  the  time  response  was  improved  in 
terms  of  overshoot  for  most  cases. 

From  these  optimizations,  a  few  important  facts  become  evident.  First,  the  closed  loop 
response  of  the  system  can  be  made  more  uniform  across  the  operating  envelope.  Second,  the 
most  effective  gain  schedule  optimization  results  when  the  control  parameters  are  scheduled 
as  piecewise  linear  functions  of  the  scheduling  variable.  Additionally,  the  number  and  size  of 
the  scheduling  variable  intervals  are  allowed  to  vary  as  design  variables.  Furthermore,  when 
using  a  central  controller  for  computation  of  gain  scheduling  error,  the  selection  of  which 
controller  is  central  can  be  an  effective  design  variable  in  reducing  the  calculated  error. 

To  see  how  the  central  controller  varied  with  the  optimization  cases  see  Fig.  5.75.  The 
‘x’  marks  the  central  controller  selected  for  the  baseline  design  using  engineering  judgment. 
The  ‘o’s’  mark  the  central  controllers  selected  by  the  optimization.  The  next  section  uses  the 
best  optimized  schedule.  Case  VEH,  and  validates  the  schedule  using  operating  points  not  in 
the  original  design  set. 

5.4.12  Design  Validation.  To  validate  the  optimal  gain  schedule  designed  in  the 
previous  section,  six  more  flight  conditions  are  selected  across  the  flight  envelope.  These  are 
shown  with  ‘o’s’  in  Fig.  5.76.  These  operating  points  are  specifically  selected  to  be  near 
the  comers  of  the  operating  envelope.  The  relative  error  and  time  response  of  these  points 
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Table  5.2  Comparison  of  Optimization  Results 
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Figure  5.75  Central  Controller  Variations 
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Figure  5.76  F- 1 8  Flight  Envelope 


using  the  baseline  schedule  are  shown  in  Fig.  5.77  and  Fig.  5.78.  The  relative  error  and  time 
response  of  these  points  using  the  optimized  schedule  are  shown  in  Fig.  5.79  and  Fig.  5.80. 

The  purpose  of  this  step  is  to  complete  the  gain  scheduling  design  process  for  the 
optimally  selected  schedules.  By  inspecting  the  relative  error  graphs,  Fig.  5.77  and  Fig. 
5.79,  all  the  relative  errors  are  less  than  one  therefore  stability  is  maintained.  The  difference 
is  that  for  the  baseline  schedule  the  relative  error  for  these  six  flight  conditions  is  2.87, 
whereas,  for  the  optimal  schedule  the  relative  error  is  only  1.1.  For  both  schedules,  the 
time  responses  have  acceptable  characteristics.  However,  with  a  lower  relative  error  of 
the  optimally  scheduled  controller,  it  is  easier  to  design  an  outer-loop  controller  to  provide 
robust  performance  throughout  the  operating  envelope.  From  the  lessons  learned  in  these 
optimizations,  a  general  design  procedure  is  developed  in  the  next  chapter. 
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Figure  5.77  Baseline  Relative  Error  Validation 
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Figure  5.79  Optimial  Schedule  Relative  Error  Validation 


Figure  5.80  Optimal  Schedule  Time  Response  Validation 
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VI.  General  Design  Process 


A  great  deal  of  experience  was  gained  in  designing  and  optimizing  the  gain  schedules 
presented  in  chapters  IV  and  V.  From  this  experience  a  systematic  gain  scheduling  method 
has  been  developed.  The  next  section  presents  this  formal  design  process. 

6.1  Gain  Scheduling  Design 

There  are  many  aspects  of  the  system  to  consider  when  designing  a  controller  for  a 
real  world  system.  Traditionally,  the  method  of  implementing  the  designed  controller  is  not 
considered  until  the  final  stages  of  the  design  process.  Ideally,  the  designed  controller  should 
be  robust  enough  to  not  require  gain  scheduling,  but  more  often  than  not  some  scheduling  is 
required.  The  design  method  presented  here  differs  fundamentally  from  most  design  methods 
in  that  the  scheduling  of  the  controller  is  considered  at  the  beginning  of  the  design  process. 

There  are  four  steps  in  the  design  process  and  are  as  follows: 

I.  Control  Problem  Definition  First  the  control  designer  must  select  the  scheduling  variable 

based  on  criteria  provided  in  [41, 30, 43, 46].  The  system  should  be  evaluated  to  avoid 
selecting  a  poor  scheduling  variable  as  outlined  in  [43,  46].  Second,  the  family  of 
plants  is  chosen  to  represent  the  extrema  of  the  operating  envelope  and  to  sufficiently 
represent  the  areas  between  these  extrema.  The  plants  should  represent  the  plant 
dynamics  accurately.  Additionally,  for  areas  of  the  operating  envelope  where  the  plant 
dynamics  change  rapidly,  additional  plants  should  be  chosen. 

II.  Controller  Design  Next  select  a  general  form  of  the  desired  controller.  After  defining  the 

form  of  the  controller,  determine  the  number  of  controller  parameters,  nscp,  that  are  to 
be  scheduled.  For  some  controllers  all  of  the  parameters  are  scheduled,  such  as  a  full 
state  feedback  gain  matrix.  In  other  controllers,  such  as  the  one  in  Chapter  V,  only  some 
of  the  controller  parameters  vary.  Once  the  scheduled  parameters  are  identified,  design 
a  controller  at  the  extrema  of  the  scheduling  variable  to  meet  the  requirements  of  the 


6-1 


system.  This  provides  upper  and  lower  bound  of  the  scheduled  variables,  although  the 
control  designer  can  extend  these  bounds  to  increase  the  search  space.  By  bounding  the 
parameters  the  control  designer  is  able  to  define  the  space  of  the  optimal  gain  schedule 
controller. 

III.  Gain  Schedule  Optimization  Next  the  gain  schedule  is  optimized  using  GAs.  The 

scheduling  variable  is  decomposed  into  n  intervals  with  each  interval  having  an  arbitrary 
size.  The  scheduled  parameters  are  bound  by  the  designs  performed  in  Step  11  and 
allowed  to  vary  linearly  in  each  interval.  The  number  of  intervals  is  chosen  to  vary  within 
a  selected  range,  n  G  {nm.in,nraax)-  Additionally,  the  central  flight  condition  is  allowed 
to  vary  within  the  family  of  plants  chosen  in  Step  I.  For  the  optimization,  the  number  of 
design  variables  varies  between  {rimin  +  x  and{nmax  +  ^max  x  Uscp).  An 
objective  function  is  then  defined  that  measures  the  gain  scheduling  error  of  the  family 
of  plants.  Now,  a  GA  can  be  used  to  optimize  the  objective  function  with  the  design 
variables  described. 

IV.  Analyze  Results  Finally,  select  a  family  of  plants  not  included  in  the  family  of  Step  I. 

Simulate  the  closed  loop  response  of  the  system  with  the  scheduled  controller  to  insure 
that  the  system  requirements  are  achieved. 

6.2  Conclusions 

The  objectives  set  forth  in  Chapter  I  were  achieved.  First,  a  simple  GA  was  found 
that  could  be  implemented  with  a  common  computer-aided  control  design  software  package. 
This  enables  the  designer  to  quickly  and  efficiently  analyze  and  optimize  a  gain  scheduled 
controller.  Second,  a  method  of  evaluating  and  measuring  the  effectiveness  of  a  gain  schedule 
design  was  developed.  This  evaluation  measure  allowed  for  the  comparison  of  various  gain 
schedules.  Next,  using  the  gain  schedule  measurement,  a  simple  gain  scheduling  problem 
was  optimized  to  validate  the  use  of  GAs  as  an  optimization  tool.  Then,  the  gain  scheduling 
optimization  method  was  used  to  evaluate  a  gain  schedule  designed  for  an  F-18  fighter  aircraft. 
It  is  shown  in  Chapter  V  that  the  gain  schedule  originally  developed  could  be  optimized  and 
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the  resulting  time  response  of  the  system  was  improved.  Finally,  the  lessons  learned  from  this 
research  were  compiled  into  a  formal  design  method  for  designing  a  gain  scheduled  controller. 

6.3  Future  Research 

Future  research  efforts  could  focus  on  improving  the  computational  efficiency  of  the 
genetic  algorithm.  Additionally,  the  precision  of  the  GA  could  be  improved  by  using  a  floating 
point  chromosome  and  a  variable  mutation  operator.  These  improvements  focus  on  the  GA 
itself.  The  basic  approach  employed  in  this  effort  is  to  decompose  the  gain  schedule  into  a 
piecewise  function  of  the  scheduling  variable.  Future  work  could  analyze  the  direct  implement 
ion  and  optimization  of  a  polynomial  scheduling  function.  Further  research  is  also  needed  in 
the  area  of  determining  global  stability  and  performance  for  a  scheduled  controller. 
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Appendix  A.  F-18  Design  Flight  Conditions 


Following  are  the  flight  conditions  and  their  respective  plant  matricies  that  were  used 
in  the  evaluation  of  the  gain  scheduled  controller. 


Table  A.l  Longitudinal  Flight  Conditions  For  Optimization 
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Following  are  the  short  period  longitudinal  dynamic  plant  matricies.  The  general  form 
of  the  longitudinal  short  period  dynamics  is  shown  in  Eq  (A.l). 
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Table  A.2  Longitudinal  Flight  Conditions  For  Evaluation 
Mach  Number  Altitude  (ft)  q  (psf)  a  (deg) 
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1.0 

The  notation  with  the  A  and  B  matricies  denotes  the  flight  condition.  For  example, 
Am3h26  denotes  the  A  matrix  at  Mach  0.3  and  altitude  26000  feet. 


4m3/i26 

‘^long 


Am5h40  _ 

^long 


Am6h30  _ 
■^long 


/\m4h22  _ 

■^long 


Am4h6  _ 

^long 


Am7h\4 

■^long 


■^long 


Am6hl5 

•^long 


-0.2296  0.9931 

0.02436  -0.2046 

-0.2423  0.9964 

-2.342  -0.1737 

-0.5088  0.994 

-1.131  -0.2804 

-0.4285  0.9916 

-0.7473  -0.3123 

-0.8018  0.9847 

-1.521  -0.5944 

-1.175  0.9871 

-8.458  -0.8776 

-0.8930  0.9852 

-4.1582  -0.6873 

-0.9181  0.9872 

-6.2419  -0.6920 


Rm3/i26  _ 

^long 


DmSMO 

■^long 


]Dm6h30 

^long 


Drn4/i22 

^long 


jymAhG 

^long 


r}m7h\4 

^long 


DmS/ilO 

^long 


omGhXh  _ 

^long 


-0.04034  -0.01145 
-1.73  -0.517 

-0.0416  -0.01141 
-2.595  -0.8161 

-0.09277  -0.01787 
-6.573  -1.525 

-0.0813  -0.0145 
-4.0770  -0.7978 

-0.1508  -0.02776 
-7.926  -1.751 

-0.194  -0.04349 
-19.29  -3.803 

-0.1626  -0.0261 
-10.6454  -2.0318 

-0.1599  -0.0265 
-12.5295  -2.4169 


A-2 


AmlhlS.h  _ 

■^long 


-0.9920 

0.9888 

jDmThlS.S  _ 

^long  “* 

-0.1652 

-0.0274 

-7.8450 

-0.7525 

-16.1838 

-2.8619 

AmQh2  _ 

^long 


-1.4710 

0.9808 

dto6/i2  _ 

^long 

-0.2419 

-0.0417 

-11.5022 

-1.0846 

-21.5227 

-3.9738 

AmShlA 

'^long 


-1.4406  0.9868 

-14.2709  -1.0645 


r>m8hl4 

^long 


-0.2154  -0.0379 
-24.3921  -4.5141 


Am8hl2 

^long 


-1.562  0.9862 

-14.94  -1.132 


r>m8h\2  _ 

^long 


-.2316  -.04349 
-26.48  -5.323 


Am95h20  _ 
^long 


-1.905 

0.9895 

E>m95/i20  _ 

^long 

-0.1867 

-0.03287 

-33.88 

-0.9872 

-27.22 

-4.573 

Am&hW  _ 
■^long 


-1.675  0.9853 

-16.16  -1.212 


jDmShlO  _ 

■^long 


-0.2449  -0.04649 
-28.34  -5.742 


jm9/il4  _ 

^long 


-2.1163 

-32.6459 


0.9872 

-1.1826 


dto9/i14 

^long 


-0.2450  -0.0426 
-32.6358  -5.7862 


AmShS  _ 
■^long 


-1.994  0.9828 

-19.44  -1.427 


Dm8/i5  _ 

^long 


-0.2852  -0.05567 
-33.44  -6.931 


/1?7i9/i10  _ 

'^long 


-2.452  0.9856 
-38.61  -1.34 


Dm9/il0 

^long 


-0.2757  -0.05226 
-37.36  -7.247 


A  7Tl85/l5 
-^long 


-2.328  0.9831 

-30.44  -1.493 


iDmSShS  _ 

^long 


-0.3012  -0.05866 
-38.43  -7.815 


^long 


-2.8375  0.9855 

-51.8325  -1.4037 


TDni95h9 

^long 


-0.2863  -0.0454 
-42.9285  -6.6039 


Am9h5 

•^long 


-2.911 

0.9835 

O  Ai  ^ 

-0.3161 

-0.06231 

i  y  'iLjilaj 

^long 

-46.47 

-1.553 

-43.65 

-8.752 

The  following  flight  conditions  are  those  used  to  evaluate  the  optimal  gain  schedule. 
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4m98/i40 

'^long 


Am99hl0 

'^long 


Am5h20  _ 

^long 


AmZhib 

■^long 


Am,2hl  _ 

^long 


AmSh\  _ 
^long 


-0.7055 

0.9949 

-17.5821 

-0.4583 

-2.6317 

0.9860 

-76.1833 

-1.3868 

-0.6199 

0.9900 

-1.8909  - 

-0.4433 

-0.3664 

0.9887 

0.0321  - 

-0.3425 

-0.3678  0.9843 

0.3351  -0.3022 

-2.2679  0.9803 

-22.8644  -1.6265 


D?n98fe40  _ 

^long 


Dm3hl5 

^long 


r)m2hl  _ 

^long 


DmShl  _ 

^long 
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Appendix  B.  Computer  Codes 


B.l  Example  Problem  Code 

Following  is  the  Matlab®  code  that  was  used  to  evaluate  the  objective  function  for 
the  sample  control  problem  in  chapter  IV.  There  are  four  m-files  shown.  The  are  the  fixed 
interval  objective  function,  the  variable  interval  size  objective  function,  the  variable  interval 
number  and  size  objective  function,  and  the  linear  schedule  objective  function,  respectively. 

%******************************************  file  fixint.m  **** 
function  [f ]  =fixint  (k)  ; 

%  This  function  is  an  evaluation  routine  run  from  genesis  to 
%  evaluate  the  distance  from  the  closed  loop  poles  of  a  system 
%  with  a  given  gain  to  the  desired  closed  loop  poles  over 
%  a  range  of  a  varying  parameter  c. 

%  The  variable  k  is  a  vector  of  gains  for  each  interval. 

%  Define  plant  transfer  function  in  form  1/  (s+a)  (s"2  +  bs  +  c) 

%  Controller  k  is  proportional 

a=6;b=4; 

%  The  characteristic  equation  of  the  closed  loop  system  is 
%  s''3  +  (a+  b)  s''2  +  (ab+c)  s  +  ac  +  k 

total=0; 

interval=[0  2468  10]; 
for  i=l:5 

for  c  =  interval (i) : 0 . 02 : interval (i+1) 
r  =  roots ([1,  (a+b),  (a*b+c) ,  (a*c+k (i) ) ] ) ' ; 
temp  =  (max  (imag  (r)  )  -  2)  "2  +  (max  (real  (r)  )  +  2)  "2; 
total  =  total  +  temp; 
end 
end 

f=total; 

return 
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^** **************************************** *  0nd.  of  file  **** 


%*****************************************  file  varint.m  **** 
function  [f ] =varint (k) ; 

%  This  function  is  an  evaluation  routine  run  from  genesis  to 
%  evaluate  the  distance  from  the  closed  loop  poles  of  a 
%  system  with  a  given  gain  to  the  desired  closed  loop 
%  poles  over  a  range  of  a  varying  parameter  c. 

%  The  variable  k  is  a  vector  of  gains  for  each  interval 
%  and  the  four  interval  break  points. 

%  Define  plant  transfer  function  in  form  1/  (s+a)  (s''2  +  bs  +  c) 
%  Controller  k  is  proportional 

a=6;b=4; 

%  The  characteristic  equation  of  the  closed  loop  system  is 
%  s'S  +  (a+  b)  s"2  +  (ab+c)  s  +  ac  +  k 

total=0; 

%  sort  the  interval  break  points  k(l)-k(4)  in  order  from 
%  lowest  to  highest 

for  1=1:3 
for  i=  1 : 3 
if  k(i)>k(i+l) 
temp  =  k (i) ; 
k(i)  =  k(i+l) ; 
k (i+1) =temp; 
end 
end 
end 

interval=[0  k(l;4)  10]; 
for  i=l:5 

for  c  =  interval (i) : 0 . 02 : interval (i+1) 
r=roots([l,  (a+b) ,  (a*b+c) ,  (a*c+k (4+i) ) ] ) ' ; 
temp  =  (max (imag (r) )  -  2) "2  +  (max (real (r) )  +  2) "2; 
total  =  total  +  temp; 
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end 

end 

f=total; 

return 

%*******************************************  end  of  file  **** 

%****************************************  file  varintn.m  **** 
function  [f ] =varintn (k) ; 

%  This  function  is  an  evaluation  routine  run  from  genesis 
%  to  evaluate  the  distance  from  the  closed  loop  poles  of 
%  a  system  with  a  given  gain  to  the  desired  closed  loop 
%  poles  over  a  range  of  a  varying  parameter  c. 

%  The  variable  k  is  a  vector  of  gains  for  each  interval 
%  and  the  four  interval  break  points.  The  first  element 
%  of  k  is  the  number  of  intervals  (2:8),  the  next  eight 
%  elements  are  the  break  points  (0:10),  and  the  final 
%  nine  elements  are  the  gains  for  each  interval  (-50:50) 

%  Define  plant  transfer  function  form  l/(s+a)  (s’'2  +  bs  +  c) 

%  Controller  k  is  proportional 

a=6;b=4; 

%  The  characteristic  equation  of  the  closed  loop  system  is 
%  s"3  +  (a+  b)  s"2  +  (ab+c)  s  +  ac  +  k 

total=0; 

%  sort  the  interval  break  points  k(l)-k(4)  in  order  from 
%  lowest  to  highest 

int=[0  k(2:k(l) )  10] ; 
s=length (int) ; 

for  1=1:3 
for  i=  1:3 
if  k(i)>k(i+l) 
temp  =  k (i) ; 
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k(i)  =  k(i+l) ; 
k (i+1) =temp; 
end 
end 
end 

interval= [0  k (1 : 4)  10]; 
for  i=l :k (1) 

for  c  =  interval (i) : 0 . 02 : interval (i+1) 
r  =  roots  ([1,  (a+b) ,  (a*b+c) ,  (a*c+k ( 9+i) ) ] ) ' ; 
temp  =  (max (imag(r) )  -  2) '2  +  (max (real  (r)  )  +  2) "2; 
total  =  total  +  temp; 
end 
end 

f=total; 

return 

%*******************************************  end  of  file  **** 

%***************************************  file  sampcoef .m  **** 
function  [f ] =sampcoef (k) ; 

%  This  function  is  an  evaluation  routine  run  from  genesis 
%  to  evaluate  the  distance  from  the  closed  loop  poles  of 
%  a  system  with  a  given  gain  to  the  desired  closed  loop 
%  poles  over  a  range  of  a  varying  parameter  c. 

%  The  variable  k  is  a  vector  of  the  coefficients  of  a 
%  polynomial  function  of  the  scheduling  variable  c. 

%  Define  plant  transfer  function  form  l/(s+a)  (s''2  +  bs  +  c) 

%  Controller  k  is  proportional 

a=6;b=4; 

%  The  characteristic  equation  of  the  closed  loop  system  is 
%  s"3  +  (a+  b)  s"2  +  (ab+c)  s  +  ac  +  k 

total=0; 
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for  c  =  0:0.02:10 
k=polyval (coef , c)  ; 

r  =  roots ([1,  (a+b) ,  (a*b+c) ,  (a*c+k)])'; 
temp  =  (max (imag (r) )  -  2)  “2  +  (max(real  (r) )  +  2) ''2; 
total  =  total  +  temp; 
end 
end 

f=total; 

return 

%*******************************************  end  of  file  **** 

B.2  F-18  Example  Code 

The  Matlah®  code  for  the  F-18  example  varies  for  each  optimization  example.  First, 
the  ‘C’  code  that  calls  the  Mat  lab©  engine  is  presented.  Next  the  code  for  evaluating  the 
baseline  design  is  presented.  Next,  for  simplicity  the  code  for  only  Cases  I,  II,  V  and  VII  are 
presented. 

%*****************************************  file  relerr.c  **** 

#include  "extern/ include/engine. h" 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <string.h> 

#include  "global. h" 

double  eval(str,  length,  vect,  genes) 

char  str[];  /*  string  representation  */ 

int  length;  /*  length  of  bit  string  */ 

double  vect[];  /*  floating  point  representation  */ 

int  genes;  /*  number  of  elements  in  vect  */ 

{ 

Matrix  *v,  *d; 
register  int  i; 
double  sum,  *Dreal; 
sum  =  0.0; 

/*  print f ("here") ; 
print  f ( " %g\n" , *vect ) ;  * / 
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V  =  mxCreateFull (1, genes, REAL) ; 

memcpy (mxGetPr (v) ,vect, genes*sizeof (double) ) ; 

mxSetName (v, "V") ; 

engPutMatrix (epl,  v) ; 
engEvalString(epl, "d=relerr (V) ") ; 
d  =  engGetMatrix (epl, "d") ; 


Dreal  =  mxGetPr(d)  ; 

/*  printf ("%g\n", *Dreal) ;  */ 

mxFreeMatrix (v)  ; 
mxFreeMatrix (d) ; 

return  (*Dreal) ; 


} 

%*******************************************  end  of  file  **** 

%*****************************************  file  relbase.m  **** 

%  This  program  calculates  the  open  and  closed  loop  dynamics  of  the  F-18 
%  short  period  longitudinal  dynamics  and  the  Relative  error  of  the  closed 
%  inner  loop  for  the  chosen  central  flight  condition. 

clear  numtemp  dentemp 
clear  nump  denp 

F  =  -40;  Kf  =  [1  1] ;  G  =  .0247; 

%  Vector  of  dynamic  pressure  for  20  flight  conditions 
%  (Original  12+8  more  off  design) 

q=[47.4  68.5  100.1  158.4  189.9  255.0  301.1  355.0  426.4  496.0  557.0  ... 
603.0  614.4  652.0  705.0  789.1  825.2  890.8  956.0  998.7] 

%  Calculate  the  scheduled  parameters  of  the  inner  loop  controller 


N=-.312*q+461 

Ml=-.058*q+50.5 

M2=-.006*q+8.11 
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%  Calculate  the  inner  loop  controller  state  space  form 


1=1; 

for  i=l : 20 
Akeq(i)=F-G*N(i) ; 

Bkeq(i, :)=Kf-G*[Ml(i)  M2(i)]; 

Ckeq(i) =-N (i) ; 

Dkeq(i, : ) =- [Ml (i)  M2(i)]; 
end 

%  The  central  controller  is  chosen  to  be  the  Am95h20 

Akcen=Akeq(13) ; 

Bkcen=Bkeq(13, :) ; 

Ckcen=Ckeq(13) ; 

Dkcen=Dkeq(13, :) ; 

%  Input  the  dynamics  at  the  20  flight  conditons 

Am3h26=[-.2296  ,9931;. 02436  -.2046]; 

Am5h40=[-.2423  .9964;-2.342  -.1737]; 

Am4h22=[-.4285  . 9916; -. 7473  -, 3123] ;  %Am4h22  *3 

Am6h30=[-.5088  .994;-!. 131  -.2804]; 

Am4h6=[-.8018  ,9847;-!. 521  -.5944]; 

Am5hl0=[-.8930  . 9852;-4 . 1582  -. 6873] ;  %Am5hl0  *6 
Am6hl5=[-.9181  . 9872;-6. 2419  -. 6920] ;  %Am6hl5  *7 
Am7hl8_5=[-.9920  .  9888; -7 . 8450  - .7525] ;  %Am7hl8_5  *8 
Am7hl4= [-1.175  .9871;-8.458  -.8776]; 

Am6h2= [-1.4710  . 9808; -11 . 5022  -1 . 0846] ;  %Am6h2  *10 
Am8hl4=[-2.1163  . 9872; -32 . 6459  -1 . 1826] ;  %Am8hl4  *11 
Am8hl2= [-1.562  .9862;-14.94  -1.132]; 

Am95h20= [-1,905  .9895;-33.88  -.9872]; 

Am8hl0= [-1.675  .9853;-16.16  -1.212]; 

Am9hl4=[-2.1163  ,  9872;-32 . 6459  -1.1826] ;  %Am9hl4  *15 
Am8h5=[-1.994  .9828;-19.44  -1.427]; 

Am9hl0=[-2.452  .9856;-38.61  -1.34]; 

Am95h9=[-2.8375  . 9855; -51 . 8325  -1 . 4037] ;  %Am95h9  *19 

Am85h5=[-2.328  .9831;-30.44  -1.493]; 

Am9h5=[-2.911  .9835;-46.47  -1.553]; 

Acen=[-1.905  .9895;-33.88  -.9872]; 

C=[l  0;0  1]; 


B-7 


D=[0;0]; 

Blong=[0;l] ; 

%  The  elevator  actuator  dynamics 
numelev= [1/82 . 9"2  2*. 068/82. 9  1]; 

denelev=conv{ [1/36.4^2  2*. 41/36. 4  1], [1/105.3^2  2*. 59/105. 3  1]) ; 

[ aact ,  bact , cact , dact ] =t  f 2ss (numelev, denelev) ; 

%break 

%  Combine  the  actuator  dynamics  in  the  plant 

[Apia, Bpla, Cpla, Dpla] =series (aact , bact , cact , dact , Am3h2  6, Blong, C, D) ; 
[Ap2a,Bp2a, Cp2a,Dp2a] =series (aact, bact, cact, dact, Am5h40, Blong, C,D) ; 
[Ap3a,Bp3a, Cp3a,Dp3a] =series ( aact, bact, cact, dact, Am4h22, Blong, C,D) ; 
[Ap4a,Bp4a, Cp4a,Dp4a] =series (aact, bact, cact, dact, Am6h30, Blong, C,D)  ; 
[Ap5a,Bp5a, Cp5a,Dp5a] =series ( aact, bact, cact, dact, Am4h6, Blong, C,D) ; 
[Ap6a,Bp6a, Cp6a,Dp6a] =series ( aact, bact, cact , dact, Am5hl0, Blong, C,D) ; 
[Ap7a,Bp7a, Cp7a,Dp7a] =series (aact, bact, cact, dact, Am6hl 5, Blong, C,D) ; 
[Ap8a,Bp8a, Cp8a,Dp8a] =series ( aact, bact, cact, dact, Am7hl8_5, Blong, C,D) ; 
[Ap9a,Bp9a, Cp9a,Dp9a] =series (aact,bact, cact, dact, Am7hl4, Blong, C,D) ; 
[Apl0a,Bpl0a, Cpl0a,Dpl0a] =series (aact, bact, cact, dact, Am6h2, Blong, C,D) ; 
[Aplla,Bplla, Cplla,Dplla] =series (aact, bact, cact, dact, Am8hl4, Blong, C,D) ; 
[Apl2a,Bpl2a, Cpl2a,Dpl2a] =series (aact, bact, cact, dact, Am8hl2, Blong, C,D) ; 
[Apl3a,Bpl3a, Cpl3a,Dpl3a] =series (aact, bact, cact, dact, Am95h20, Blong, C,D) ; 
[Apl4a,Bpl4a,Cpl4a,Dpl4a]=series (aact, bact, cact, dact, Am8hl0, Blong, C,D) ; 
[Apl5a,Bpl5a, Cpl5a,Dpl5a] =series (aact, bact, cact, dact, Am9hl4, Blong, C,D) ; 
[Apl 6a, Bpl 6a, Cpl 6a, Dpi 6a] =series (aact , bact , cact , dact , Am8h5, Blong, C,D); 
[Apl7a,Bpl7a, Cpl 7a, Dpi 7a] =series (aact, bact, cact, dact, Am9hl0, Blong, C,D) ; 
[Apl 8a, Bpl 8a, Cpl 8a, Dpi 8a] =series (aact, bact, cact, dact, Am85h5, Blong, C,D) ; 
[Apl 9a, Bpl 9a, Cpl9a,Dpl9a] =series (aact,bact, cact, dact, Am95h9, Blong, C,D) ; 
[Ap20a,Bp20a, Cp20a,Dp20a]=series (aact, bact, cact, dact, Am9h5, Blong, C,D) ; 
[Ap0a,Bp0a, Cp0a,Dp0a] =series (aact, bact, cact, dact, Acen, Blong, C,D) ; 

%  Close  the  inner  loop  with  the  corresponding  controller 

[ Apl, Bpl, Cpl, Dpi ]=feedback( Apia, Bpla, Cpla, Dpla, . . . 

Akeq(l) ,Bkeq(l, : ) , Ckeq(l) ,Dkeq(l,  :),1); 
[Ap2,Bp2,Cp2,Dp2]=feedback(Ap2a,Bp2a, Cp2a,Dp2a, . . . 

Akeq(2) ,Bkeq(2, : ) , Ckeq(2) ,Dkeq(2,  : ) ,  1) ; 
[Ap3,Bp3,Cp3,Dp3]=feedback(Ap3a,Bp3a, Cp3a,Dp3a,  . . . 

Akeq(3) ,Bkeq(3, : ) , Ckeq(3) ,Dkeq(3, : ) , 1) ; 
[Ap4,Bp4,Cp4,Dp4]=feedback(Ap4a,Bp4a, Cp4a,Dp4a, . . . 
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Akeq(4)  ,Bkeq(4, : ) , Ckeq(4) ,Dkeq(4, : ) , 1) ; 
[Ap5,Bp5,Cp5,Dp5]=feedback(Ap5a,Bp5a,Cp5a,Dp5a, , . . 

Akeq(5) ,Bkeq(5, :) ,Ckeq(5) ,Dkeq(5,  :)  ,1)  ; 
[Ap6,Bp6,Cp6,Dp6]=feedback(Ap6a,Bp6a, Cp6a,Dp6a, . . . 

Akeq(6) ,Bkeq(6, :) ,Ckeq(6) ,Dkeq(6,  :)  ,1) ; 
[Ap7,Bp7,Cp7,Dp7]=feedback(Ap7a,Bp7a, Cp7a,Dp7a, . . . 

Akeq(7) ,Bkeq(7, :) ,Ckeq(7)  ,Dkeq(7,  : ) ,  1)  ; 
[Ap8,Bp8,Cp8,Dp8]=feedback(Ap8a,Bp8a, Cp8a,Dp8a, , , . 

Akeq(8) ,Bkeq(8, : ) ,Ckeq(8)  ,Dkeq(8,  : ) ,  1)  ; 
[Ap9,Bp9,Cp9,Dp9]=feedback(Ap9a,Bp9a, Cp9a,Dp9a, . . . 

Akeq(9) ,Bkeq(9, :) ,Ckeq(9) ,Dkeq(9,  :)  ,1) ; 

[AplO,BplO,CplO,DplO]=feedback(AplOa,BplOa,CplOa,DplOa, , . 
Akeq(lO) ,Bkeq(10, : ) , Ckeq(lO)  ,Dkeq(10,  : ) ,  1)  ; 
[Apll,Bpll,Cpll,Dpll]=feedback(Aplla,Bplla,Cplla,Dplla, . . 
Akeq(ll) ,Bkeq(ll, : ) , Ckeq(ll) ,Dkeq(ll,  : ) ,  1) ; 
[Apl2,Bpl2,Cpl2,Dpl2]=feedback(Apl2a,Bpl2a,Cpl2a,Dpl2a, , . 
Akeq(12) ,Bkeq(12, : ) , Ckeq(12) ,Dkeq(12,  : ) ,  1) ; 
[Apl3,Bpl3,Cpl3,Dpl3]=feedback(Apl3a,Bpl3a,Cpl3a,Dpl3a, . . 
Akeq(13) ,Bkeq(13, : ) , Ckeq(13) ,Dkeq(13,  : ) ,  1) ; 

[Apl4,Bpl4, Cpl4,Dpl4] =feedback (Apl4a,Bpl4a, Cpl4a,Dpl4a,  . . 
Akeq(14),Bkeq(14,:),Ckeq(14),Dkeq(14,:),l); 
[Apl5,Bpl5,Cpl5,Dpl5]=feedback{Apl5a,Bpl5a,Cpl5a,Dpl5a, . . 
Akeq(15) ,Bkeq(15, : ) , Ckeq(15) ,Dkeq(15, 

[Apl6,Bpl6, Cpl6,Dpl6] =feedback (Apl6a,Bpl6a, Cpl6a,Dpl6a, . . 
Akeq(16)  ,Bkeq(16,  : ) , Ckeq(16)  ,Dkeq(l6, 

[Apl7,Bpl7,Cpl7,Dpl7]=feedback(Apl7a,Bpl7a,Cpl7a,Dpl7a, . . 
Akeq(17) ,Bkeq(17, : ) , Ckeq(17) ,Dkeq(17,  : ) , 1) ; 
[Apl8,Bpl8,Cpl8,Dpl8]=feedback(Apl8a,Bpl8a,Cpl8a,Dpl8a, . . 
Akeq(18) ,Bkeq(18, : ) , Ckeq(18) ,Dkeq(18, 

[Apl9,Bpl9, Cpl9,Dpl9] =feedback (Apl9a,Bpl9a, Cpl9a,Dpl9a, . . 
Akeq(19) ,Bkeq(19, : ) , Ckeq(19) ,Dkeq(19, :),1) ; 
[Ap20,Bp20,Cp20,Dp20]=feedback(Ap20a,Bp20a,Cp20a,Dp20a, . . 
Akeq(20) ,Bkeq(20, : ) , Ckeq(20) ,Dkeq(20,  : ) , 1) ; 
[ApO,BpO,CpO,DpO]=feec[back(ApOa,BpOa, CpOa,DpOa,  .  .  . 
Akcen,Bkcen, Ckcen,Dkcen,  1)  ; 

[Apl,Bpl,Cpl,Dpl]=feedback(Am3h26,Blong,C,D, . . . 

Akeq ( 1 ) , Bkeq ( 1 , : ) , Ckeq ( 1 ) , Dkeq ( 1 ,  : ) ,  1 ) ; 

[Ap2,Bp2, Cp2,Dp2] =feedback (Am5h40,Blong, C,D, . . . 

Akeq(2) ,Bkeq(2, :) ,Ckeq(2) ,Dkeq(2,  :) ,1) ; 
[Ap3,Bp3,Cp3,Dp3]=feedback(Am4h22,Blong,C,D, ... 

Akeq (3) , Bkeq (3, : ) , Ckeq (3) , Dkeq (3, : ) , 1) ; 
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[Ap4,Bp4,Cp4,Dp4]=feedback(Am6h30,Blong, C,D, . . . 

Akeq(4) ,Bkeq(4, : ) , Ckeq(4) ,Dkeq(4, 

[Ap5,Bp5, Cp5,Dp5] =feedback (Am4h6,Blong, C,D, . , . 

Akeq(5) ,Bkeq(5, :) ,Ckeq(5) ,Dkeq(5,  :)  ,1) ; 
[Ap6,Bp6,Cp6,Dp6]=feedback(Am5hlO,Blong,C,D, , . . 

Akeq(6) ,Bkeq(6, :) ,Ckeq(6) ,Dkeq(6, 
[Ap7,Bp7,Cp7,Dp7]=feedback(Am6hl5,Blong,C,D, . . . 

Akeq ( 7 ) , Bkeq ( 7 , : ) , Ckeq ( 7 ) , Dkeq (7,:),1); 

[Ap8,Bp8, Cp8,Dp8] =feedback (Am7hl8_5,Blong, C,D, . . . 

Akeq(8) , Bkeq (8, : ) , Ckeq (8)  , Dkeq (8, 

[Ap9,Bp9,  Cp9,Dp9]  =feedback  (Ain7hl4,Blong,C,D,  .  .  . 

Akeq(9) ,Bkeq(9, :) ,Ckeq(9) ,Dkeq(9, 

[AplO,BplO, CplO,DplO] =feedback (Am6h2,Blong, C,D, . , . 

Akeq(lO) ,Bkeq(10, :) ,Ckeq(10) ,Dkeq(10, : ) , 1) ; 

[Apll,Bpll, Cpll,Dpll] =feedback (Am8hl4,Blong, C,D, . . . 

Akeq (11) , Bkeq (11,  : ) , Ckeq (11) , Dkeq (11,  : ) , 1) ; 

[Apl2,Bpl2, Cpl2,Dpl2] =feedback (Am8hl2,Blong, C,D,  . . . 

Akeq(12) ,Bkeq(12, ;) ,Ckeq(12) ,Dkeq(12,  : ) , 1) ; 
[Apl3,Bpl3,Cpl3,Dpl3]=feedback(Am95h20,Blong,C,D,  .  , . 

Akeq (13) , Bkeq (13, : ) , Ckeq (13) , Dkeq (13,  : ) ,  1)  ; 
[Apl4,Bpl4,Cpl4,Dpl4]=feedback(Ara8hlO,Blong,C,D, . . . 

Akeq(14) ,Bkeq(14, :),Ckeq(14) ,Dkeq(14, : ) , 1) ; 
[Apl5,Bpl5,Cpl5,Dpl5]=feedback(Ain9hl4,Blong,C,D,  . .  . 

Akeq (15) , Bkeq (15, : ) ,  Ckeq (15) , Dkeq (15, : ) , 1) ; 

[Apl6,Bpl6, Cpl6,Dpl6] =feedback (Am8h5,Blong, C,D, . . . 

Akeq (16) , Bkeq (16, : ) , Ckeq (16) , Dkeq (16, : ) , 1) ; 

[Apl7,Bpl7, Cpl7,Dpl7] =feedback (Am9hlO,Blong, C,D,  . , . 

Akeq(17) ,Bkeq(17, :) ,Ckeq(17) ,Dkeq(17,  :),1); 

[Apl8,Bpl8, Cpl8,Dpl8] =feedback (Am85h5,Blong, C,D, . . . 

Akeq (18) , Bkeq (18, : ) , Ckeq (18) , Dkeq (18, : ) , 1) ; 

[Apl9,Bpl9, Cpl9,Dpl9] =feedback (Am95h9,Blong, C,D,  . . . 

Akeq(19) ,Bkeq(19, :) ,Ckeq(19) ,Dkeq(19, : ) , 1) ; 
[Ap20,Bp20,Cp20,Dp20]=feedback(Am9h5,Blong,C,D, . , . 

Akeq (20) , Bkeq (20, : ) , Ckeq (20) , Dkeq (20,  : ) ,  1)  ; 
[ApO,BpO,CpO,DpO]=feedback(Acen,Blong, C,D, . . . 

Akcen,Bkcen, Ckcen,Dkcen, 1) ; 

w=logspace (-2,2,200) ; 

C1=C(1, :) ;D1=0; 

C2=C(2,  :)  ;D2=0; 

%  Compute  open  loop  dynamics  Pitch  acceleration  to  angle  of  attack 
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[svlol] =sigma (Am3h26,Blong, Cl,Dl,w) ; 

[sv2ol]  =sigma (i\in5h40,Blong, Cl, Dl, w) ; 
[sv3ol]=sigma(M4h22,Blong, Cl,Dl,w) ; 
[sv4ol]=sigma(Am6h30,Blong, C1,D1, w) ; 

[svSol]  =sigina (Am4h6, Blong, Cl, Dl, w)  ; 

[sv6ol]  =sigma(Mi5hlO,Blong, C1,D1, w)  ; 

[sv7ol] =sigma (Am6hl5, Blong, C1,D1, w) ; 

[svSol] =sigma (Am7hl8_5, Blong, Cl,Dl,w) ; 

[sv9ol]  =sigma (M7hl4, Blong, C1,D1, w)  ; 

[svlOol]  =sigma  (Ain6h2, Blong,  Cl,Dl,w) ; 

[svllol] =sigma (Am8hl4,Blong,Cl,Dl,w) ; 

[svl2ol] =sigma (Am8hl2, Blong, C1,D1, w) ; 

[svl3ol]=sigma (Am95h20, Blong, Cl, Dl,w) ; 

[svl4ol]=sigma (Am8hl0, Blong, Cl, Dl,w) ; 

[svlSol]  =sigina (Am9hl4, Blong, C1,D1, w) ; 

[svl6ol] =sigma (Am8h5, Blong, C1,D1,  w) ; 

[svl7ol]  =sigrtia  (Mi9hl0, Blong,  C1,D1,  w)  ; 

[svl8ol] =sigma (Am85h5,Blong,Cl,Dl, w) ; 

[svl9ol]  =sigma (Ain95h9, Blong, Cl, Dl,w)  ; 

[sv20ol]=sigma (Am9h5,Blong,Cl,Dl, w) ; 

%  Compute  open  loop  dynamics  Pitch  acceleration  to  pitch  rate 

[svlo2] =sigma (Am3h26, Blong, C2,D2,w) ; 

[sv2o2] =sigma (Am5h40, Blong, C2,D2,w) ; 

[sv3o2] =sigma (Am4h22, Blong, C2,D2,w) ; 

[sv4o2] =sigma (Am6h30, Blong, C2,D2,w) ; 

[sv5o2] =sigma (Am4h6, Blong, C2,D2,w) ; 

[sv6o2] =sigma (AmShlO, Blong, C2,D2,w) ; 

[sv7o2] =sigma (Am6hl5, Blong, C2,D2, w) ; 

[ sv8o2 ]=sigma(Am7hl8_5, Blong, C2,D2,w) ; 

[sv9o2] =sigma (Am7hl4, Blong, C2,D2, w) ; 

[svl0o2] =sigma (Am6h2,Blong,C2,D2, w) ; 

[svllo2] =sigma (Am8hl4,Blong,C2,D2,w) ; 

[svl2o2] =sigma (Am8hl2, Blong, C2,D2,w) ; 

[svl3o2] =sigma (Am95h20, Blong, C2,D2,w) ; 

[svl4o2]=sigma ( Am8hl0, Blong, C2,D2,w) ; 

[ svl5o2 ] =sigma (Am9hl 4 , Blong, C2 , D2 , w) ; 

[svl6o2] =sigma (Am8h5, Blong, C2,D2, w) ; 

[svl7o2] =sigma (Am9hl 0, Blong, C2,D2,w) ; 

[svl8o2] =sigma (Am85h5, Blong, C2,D2,w) ; 
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[ svl 9o2 ]  =sigma  (Ain95h9 , Blong,  C2 , D2,  w) ; 

[sv20o2] =sigma (Am9h5, Blong, C2,D2,w) ; 

%  Closed  loop  singular  value  analysis  (pitch  acceleration  to  pitch  rate) 


[svlp2]=sigma(Apl,Bpl,Cpl (2, :) ,Dpl (2,  :) 
[sv2p2] =sigma (Ap2,Bp2, Cp2 (2, : ) ,Dp2 (2, : ) 
[sv3p2]=sigma(Ap3,Bp3,Cp3 (2, :) ,Dp3(2, :) 
[sv4p2] =sigma (Ap4,Bp4, Cp4 (2, : ) ,Dp4 (2, : ) 
[sv5p2] =sigma (Ap5,Bp5, Cp5 (2, : ) ,Dp5 (2, : ) 
[sv6p2] =sigma (Ap6,Bp6, Cp6 (2, : ) ,Dp6 (2, : ) 
[sv7p2]=sigma(Ap7,Bp7,Cp7 (2, : ) ,Dp7 (2, : ) 
[sv8p2] =sigma (Ap8,Bp8, Cp8 (2, : ) ,Dp8 (2, : ) 
[sv9p2]=sigma(Ap9,Bp9,Cp9(2, :) ,Dp9(2, :) 
[svl0p2] =sigma (AplO,BplO, CplO (2, : ) ,DplO 
[svllp2] =sigma (Apll,Bpll, Cpll (2, : ) ,Dpll 
[svl2p2] =sigma (Apl2,Bpl2, Cpl2 (2, : ) ,Dpl2 
[svl3p2] =sigma (Apl3,Bpl3, Cpl3 (2, : ) ,Dpl3 
[svl4p2]=sigma(Apl4,Bpl4,Cpl4 (2, : ) ,Dpl4 
[svl5p2] =sigma (Apl5,Bpl5, Cpl5 (2, : ) ,Dpl5 
[svl6p2] =sigma (Apl6,Bpl6, Cpl6 (2, : ) ,Dpl6 
[svl7p2]=sigma(Apl7,Bpl7,Cpl7 (2, :) ,Dpl7 
[svl8p2] =sigma (Apl8,Bpl8, Cpl8 (2, : ) ,Dpl8 
[svl9p2] =sigma (Apl9,Bpl9, Cpl9 (2, : ) ,Dpl9 
[sv20p2] =sigma (Ap20,Bp20, Cp20 (2, : ) ,Dp20 

%  Reduce  the  closed  loop  system  to  SISO 

Cpl=Cpl(l, :) ;Dpl=Dpl(l,  :)  ; 

Cp2=Cp2 ( 1 , : ) ; Dp2=Dp2 ( 1 ,  : )  ; 

Cp3=Cp3  ( 1 ,  ; )  ;  Dp3=Dp3  ( 1 ,  : )  ; 

Cp4=Cp4(l, :) ;Dp4=Dp4(l, :) ; 

Cp5=Cp5 (1, : ) ;Dp5=Dp5 (1,  :)  ; 

Cp6=Cp6 (1, : ) ;Dp6=Dp6 (1,  : ) ; 

Cp7=Cp7(l, :) ;Dp7=Dp7(l,  :)  ; 

Cp8=Cp8(l, :) ;Dp8=Dp8(l,  :)  ; 

Cp9=Cp9  (1,  : )  ;Dp9=Dp9  (1,  : )  ; 

CplO=CplO(l, :) ;DplO=DplO(l, :) ; 
Cpll=Cpll(l, :) ;Dpll=Dpll(l, :) ; 
Cpl2=Cpl2(l,  :)  ;Dpl2=Dpl2(l,  :)  ; 
Cpl3=Cpl3(l, :) ;Dpl3=Dpl3(l, :) ; 
Cpl4=Cpl4(l, :) ;Dpl4=Dpl4(l,  :)  ; 


,w)  ; 

/W)  ; 

,w)  ; 

/w)  ; 

/w)  ; 

/W)  ; 

,w)  ; 

/w)  ; 

,w)  ; 

(2, :),w); 

(2,  :),w); 

(2, : ) ,w) ; 

(2, : ) ,w) ; 

{2,:),w); 

(2,:),w); 

(2, : ) ,w) ; 

(2, :),w); 

(2, : ) ,w) ; 

(2, : ) ,w) ; 

(2, : ) ,w) ; 

(pitch  acceleration  to  angle  of  attack) 
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Cpl5=Cpl5 

(1, 

) ;Dpl5=Dpl5 

(1, 

)  ; 

Cpl6=Cpl6 

(1, 

) ;Dpl6=Dpl6 

(1, 

); 

Cpl7=Cpl7 

(1, 

) ;Dpl7=Dpl7 

(1, 

); 

Cpl8=Cpl8 

(1, 

) ;Dpl8=Dpl8 

(1, 

); 

Cpl9=Cpl9 

(1, 

) ;Dpl9=Dpl9 

(1, 

); 

Cp20=Cp20 

(1, 

) ;Dp20=Dp20 

(1, 

); 

CpO=CpO(l, :) ;DpO=DpO(l,  :) ; 

%  Closed  loop  singular  value  analysis  (pitch  acceleration  to  angle  of  attack) 

[svlp] =sigma(Apl,Bpl, Cpl,Dpl, w) ; 

[sv2p] =sigma (Ap2,Bp2, Cp2,Dp2, w) ; 

[sv3p] =sigma (Ap3,Bp3, Cp3,Dp3, w) ; 

[sv4p]  =sign:ia (Ap4,Bp4, Cp4,Dp4, w)  ; 

[sv5p] =sigma (Ap5,Bp5, Cp5,Dp5,w) ; 

[sv6p] =sigma (Ap6,Bp6, Cp6,Dp6, w) ; 

[  sv7p  ]  =sigtna  ( Ap7 ,  Bp7 ,  Cp7 ,  Dp7 ,  w)  ; 

[sv8p] =sigma (Ap8,Bp8,Cp8,Dp8, w) ; 

[sv9p] =sigma (Ap9,Bp9, Cp9,Dp9, w) ; 

[svlOp] =sigma(AplO,BplO,CplO,DplO,w) ; 

[svllp] =sigma(Apll,Bpll,Cpll,Dpll, w) ; 

[svl2p] =sigma (Apl2,Bpl2, Cpl2,Dpl2, w) ; 

[svl3p] =sigma (Apl3,Bpl3,Cpl3,Dpl3,w) ; 

[svl4p]=sigina (Apl4,Bpl4,Cpl4,Dpl4,w) ; 

[svlSp] =sigma (Apl5,Bpl5,Cpl5,Dpl5,w) ; 

[svl6p]  =sigina  (Apl6,Bpl6,  Cpl6,Dpl6,  w)  ; 

[svl7p] =sigma (Apl7,Bpl7, Cpl7,Dpl7, w) ; 

[svl8p] =sigma (Apl8,Bpl8, Cpl8,Dpl8, w) ; 

[svl9p] =sigma (Apl9,Bpl9, Cpl9,Dpl9, w) ; 

[sv20p] =sigma(Ap20,Bp20,Cp20,Dp20,w) ; 

%  Convert  state  space  to  transfer  function  for  computation  of  relative  error 

[nump(l,  :)  ,denp(l,  :)  ]=ss2tf  (Apl,Bpl,Cpl,Dpl)  ; 

[nump(2, :) ,denp(2, :) ]=ss2tf (Ap2,Bp2,Cp2,Dp2) ; 

[nump(3, :) ,denp(3, :) ] =ss2tf (Ap3,Bp3, Cp3,Dp3) ; 

[nump(4, :) ,denp(4, :) ] =ss2tf (Ap4,Bp4,Cp4,Dp4) ; 

[nump(5, :) ,denp(5, :) ] =ss2tf (Ap5,Bp5,Cp5,Dp5) ; 

[nump(6,  :)  ,denp(6,  :)  ]  =ss2tf  (Ap6,Bp6,Cp6,Dp6) ; 

[nump(7, :) ,denp(7, :) ] =ss2tf (Ap7,Bp7,Cp7,Dp7) ; 

[nump (8, : ) , denp (8, : ) ] =ss2tf (Ap8,Bp8,Cp8,Dp8) ; 

[nump ( 9, : ) ,  denp (9, : ) ] =ss2tf (Ap9,Bp9,Cp9,Dp9) ; 
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[numpdO,  ;)  ,denp(10,  :)  ]  =ss2tf  (AplO,BplO, CplO,DplO)  ; 

[numpdl,  ;)  ,denpdl,  : )  ]  =ss2tf  (Apll,Bpll,  Cpll,Dpll)  ; 

[nurapd2,  ;)  ,denpd2,  :)  ]=ss2tf  (Apl2,Bpl2,Cpl2,Dpl2)  ; 

[numpdS,  :)  ,denpd3,  : )  ]=ss2tf  (Apl3,Bpl3,Cpl3,Dpl3)  ; 

[nunip(14,  :)  ,denp(14,  :)  ]  =ss2tf  (Apl4,Bpl4, Cpl4,Dpl4)  ; 

[nunip(15,  :)  ,denp(15,  :)  ]=ss2tf  (Apl5,Bpl5,Cpl5,Dpl5)  ; 

[nump (16, : ) , denp (16, : ) ] =ss2tf (Apl6,Bpl6, Cpl6,Dpl6) ; 

[nump(17, :) ,denp(17, : ) ] =ss2tf (Apl7,Bpl7, Cpl7,Dpl7) ; 

[nunip(18,  :)  ,denp(18,  ;)  ]=ss2tf  (Apl8,Bpl8,Cpl8,Dpl8) ; 

[nump (19, ; ) , denp (19, ; ) ] =ss2tf (Apl9,Bpl9, Cpl9,Dpl9) ; 

[nuiiip(20,  ;)  ,denp(20,  :)  ]=ss2tf  (Ap20,Bp20,Cp20,Dp20)  ; 

[numpO, denpO] =ss2tf (ApO,BpO, CpO,DpO) ; 

%  Reduce  the  transfer  fuction  to  essential  parts  (chop  off  added  zeros) 
nump=nump ( : , 3 : 4 ) 
numpO=numpO (3:4) 

numpO=nump (13, : ) ; denpO=denp (13,  : ) ; 

%  Compute  Relative  Error  (Po-P)Po''-l 

for  i=l:20 

numtemp(i, : ) =conv (nump (i, :) ,denpO) -conv (numpO, denp (i, :) ) ; 
dentemp (i, : ) =conv(denp (i, : ) , numpO) ; 

end 


%  Convert  back  to  state  space  for  singular  value  analysis 

[al,bl,cl, dl]=tf2ss (numtemp (1, : ) , dentemp (1, : ) ) ; 
[a2,b2,c2, d2] =tf2ss (numtemp (2, : ) , dentemp (2, : ) ) ; 
[a3,b3,c3, d3] =tf2ss (numtemp (3, : ) , dentemp (3, :)); 
[a4,b4,c4, d4] =tf2ss (numtemp (4, : ) , dentemp (4,:)); 
[a5,b5,c5,d5]=tf2ss (numtemp (5, ;) , dentemp (5,  : ) ) ; 
[a6,b6,c6, d6]=tf2ss (numtemp (6, : ) ,  dentemp (6,  : ) )  ; 
[a7,b7,c7,d7]=tf2ss(numtemp(7, :) ,dentemp(7,  :)); 
[a8,b8,c8,d8]=tf2ss(numtemp(8, :) ,dentemp(8, : ) ) ; 
[a9,b9,c9,d9]=tf2ss (numtemp (9, : ) , dentemp (9, : ) ) ; 

[alO,blO, clO, dlO] =tf2ss (numtemp (10, :) ,dentemp(10, :) ) ; 
[all, bll, cll, dll] =tf2ss (numtemp (11, :) ,dentemp(ll,  :)); 
[al2,bl2, cl2, dl2] =tf2ss (numtemp (12, :) ,dentemp(12, : ) ) ; 
[al3,bl3, cl3, dl3] =tf2ss (numtemp (13, : ) , dentemp (13, : ) ) ; 


B-14 


[al4,bl4,cl4,  dl4] =tf2ss (numtemp(15, 
[al5,bl5, cl5, dl5] =tf2ss (numtemp (15, 
[al6,bl6, cl6,dl6] =tf2ss (numtemp (16, 
[al7,bl7, cl7,dl7] =tf2ss (numtemp (17, 
[al8,bl8, cl8,  dl8] =tf2ss (numtemp (18, 
[al9,bl9, cl9,dl9] =tf2ss (numtemp (19, 
[a20,b20, c20, d20] =tf2ss (numtemp (20, 


) , dentemp (15, 

)); 

) , dentemp (15, 

)); 

) , dentemp (16, 

)); 

) , dentemp (17, 

)); 

) , dentemp (18, 

)); 

) , dentemp (19, 

)); 

) , dentemp (20, 

)); 

%  Singular  value  analysis 

[svl]  =sigma (al,bl, cl, dl, w)  ; 

[sv2]  =sigma (a2,b2, c2, d2, w) ; 

[sv3] =sigma (a3,b3, c3, d3, w) ; 

[sv4]  =sigma (a4,b4, c4, d4, w)  ; 

[sv5] =sigma (a5,b5, c5, d5, w) ; 

[sv6] =sigma (a6,b6, c6, d6, w) ; 

[sv7] =sigma (a7,b7, c7, d7,w) ; 

[sv8] =sigma (a8,b8, c8, d8, w) ; 

[sv9] =sigma (a9,b9, c9, d9,w) ; 
[svlO]=sigma (al0,bl0, clO,  dl0,w)  ; 
[svll]=sigma (all,bll, cll, dll, w) ; 
[svl2]==sigma  (al2,bl2,  cl2,  dl2,  w) ; 
[svl3] =sigma (al3,bl3,  cl3,  dl3,  w) ; 
[svl4] =sigma (al4,bl4, cl4,  dl4,w) ; 
[svl5] =sigma (al5,bl5, cl5, dl5, w) ; 
[svl6] =sigma (al6,bl6, cl6, dl6, w) ; 
[svl7] =sigma (al7,bl7, cl7, dl7, w)  ; 
[svl8] =sigma (al8,bl8, cl8, dl8, w) ; 
[svl9] =sigma (al9,bl9, cl9, dl9, w) ; 
[sv20] =sigma (a20,b20, c20, d20, w) ; 


j=max ( svl ) +max ( sv2 ) +max ( sv3 ) +max ( sv4 ) +max ( sv5 ) +max ( sv6 ) +max ( svl ) +max ( sv8 ) + 
max (sv9) +max (svlO) +max (svll) +max(svl2) +max(svl3) +max (svl4) +max (svl5) +. . . 
max ( s vl 6 ) +max ( svl 7 ) +max ( svl 8 ) +max ( svl 9 ) +max ( sv2  0 ) 
toe 

figure  ( 1 ) 

loglog(w,  svl,  w,  sv2,w,  sv3,  w,  sv4,w,  sv5,  w,  sv6,  w,  sv7,w,  sv8,  w,  sv9,  w,  svlO,  w,svll 
w,  svl2,  w,  svl3,  w,  svl4,  w,  svl5,w,  svl6,w,  svl7,  w,  svl8,w,  svl9,  w,  sv20) 

%title ('Baseline  Closed  Inner  Loop  Relative  Error') 
xlabel('w  (rad/sec)') 
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ylabel ( '  Singular  Value' ) 
axis([.01  100  .0001  1]) 
print  -deps  relbase 

figure  (2) 

loglog(w,  svlp,  w,  sv2p,  w,  sv3p,w,  sv4p,w,  sv5p,  w,  sv6p,w,  sv7p,  w,  sv8p,  w,  sv9p, .  . . 
w,  svlOp,  w,  svllp, w, svl2p, w, svl3p,w, svl4p, w, svlSp, w, svl6p, w, svl7p, w, svl8p, . 
w, svl9p, w, sv20p) 
xlabell'w  (rad/sec)') 
ylabel ( '  Singular  Value' ) 

%title ('Pitch  Acceleration  to  Angle  of  Attack  Closed  Inner-Loop  Dynamics') 
axis([.01  100  .0001  .1]) 
print  -deps  basecla 

figure  (3) 

loglog (w,  svlol, w,  sv2ol,  w,  sv3ol, w, sv4ol, w, svSol, w, sv6ol, w, sv7ol, w, sv8ol, . . 
w,  sv9ol,  w, svlOol, w, svllol, w, svl2ol, w, svl3ol , w, svl4ol , w, svl5ol, w, svl6ol, . . 
w, svl7ol, w, svl8ol, w, svl9ol, w, sv20ol) 
xlabel('w  (rad/sec)') 
ylabel  ( '  Singular  Value' ) 

%title ('Pitch  Acceleration  to  Angle  of  Attack  Open  Inner-Loop  Dynamics') 
print  -deps  olalpha 

figure  ( 4 ) 

loglog (w, svlo2,w, sv2o2, w, sv3o2,w, sv4o2,w, sv5o2, w, sv6o2, w, sv7o2, w, sv8o2, . . 
w,  sv9o2,  w, svl0o2, w, svllo2,w, svl2o2,w, svl3o2, w, svl4o2,w, svl5o2,w, svl6o2, . . 
w, svl7o2, w, svl8o2, w, svl9o2, w, sv20o2) 
xlabel('w  (rad/sec)') 
ylabel (' Singular  Value' ) 

%title('Pitch  Acceleration  to  Pitch  Rate  Open  Inner-Loop  Dynamics') 
print  -deps  olpitch 

figure  (5) 

loglog (w,  svlp2, w, sv2p2, w, sv3p2, w, sv4p2, w, sv5p2, w, sv6p2, w, sv7p2, w, sv8p2, . . 
w,  sv9p2, w, svl0p2, w, svllp2, w, svl2p2, w, svl3p2, w, svl4p2, w, svl5p2, w, svl6p2, . . 
w, svl7p2,w, svl8p2,w, svl9p2,w, sv20p2) 
xlabel('w  (rad/sec)') 
ylabel ( ' Singular  Value' ) 

%title ('Pitch  Acceleration  to  Pitch  Rate  Closed  Inner-Loop  Dynamics') 
axis([.01  100  .001  1]) 
print  -deps  baseclp 
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%  Load  the  outer  loop  controller 
longout 

%  Combine  the  outer  loop  controller  and  the  Closed  inner  loop 

[Aool, Bool, Cool, Dool] =series (Akoutlow,Bkoutlow, Ckoutlow,Dkoutlow, . . . 
Apl,Bpl, Cpl,Dpl) ; 

[ Aoo2 ,  Boo2 , Coo2 , Doo2 ] =series ( Akoutlow, Bkout low, Ckout low, Dkout low,  . . . 
Ap2,Bp2,Cp2,Dp2) ; 

[Aoo3,Boo3, Coo3,Doo3] =series (Akoutlow, Bkout low, Ckout low, Dkout low,  . .  . 
Ap3,Bp3,Cp3,Dp3) ; 

[Aoo4,Boo4, Coo4,Doo4] =series (Akoutlow, Bkout low, Ckout low, Dkout low, . . , 
Ap4,Bp4,Cp4,Dp4) ; 

[Aoo5,Boo5,Coo5,Doo5]=series(Akouthigh,Bkouthigh,Ckouthigh,Dkouthigh, , . . 

Ap5,Bp5, Cp5,Dp5) ; 

[Aoo6,Boo6, Coo6,Doo6] =series (Akouthigh,Bkouthigh, Ckouthigh,Dkouthigh, . . , 
Ap6,Bp6,Cp6,Dp6)  ; 

[Aoo7,Boo7,Coo7,Doo7]=series(Akouthigh,Bkouthigh,Ckouthigh,Dkouthigh, . . , 
Ap7,Bp7, Cp7,Dp7) ; 

[Aoo8,Boo8,Coo8,Doo8]=series(Akouthigh,Bkouthigh,Ckouthigh,Dkouthigh, . . . 
Ap8,Bp8,Cp8,Dp8) ; 

[Aoo9,Boo9,Coo9,Doo9]=series(Akouthigh,Bkouthigh,Ckouthigh,Dkouthigh, . . . 
Ap9,Bp9,Cp9,Dp9)  ; 

[AoolO,BoolO, Cool 0, Dool 0] =series (Akouthigh,Bkouthigh,Ckouthigh,Dkouthigh, 
Apl0,Bpl0,Cpl0,Dpl0) ; 

[Aooll,Booll, Cooll,Dooll] =series (Akouthigh,Bkouthigh,Ckouthigh,Dkouthigh, 
Apll,Bpll,Cpll,Dpll) ; 

[A00I2 ,  B00I2 , C00I2 , D00I2 ] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
Apl2,Bpl2,Cpl2,Dpl2) ; 

[Aool3,Bool3, Cool3,Dool3] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
Apl3,Bpl3,Cpl3,Dpl3) ; 

[Aool4,Bool4, Cool4,Dool4] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
Apl4,Bpl4,Cpl4,Dpl4) ; 

[Aool5,Bool5, Cool5,Dool5] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
Apl5,Bpl5, Cpl5,Dpl5) ; 

[Aool 6 , Bool  6 , Cool  6 , Dool 6 ] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
Apl6,Bpl6,Cpl6,Dpl6) ; 

[ A00I7 , B00I7 , C00I7 , D00I7 ] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
Apl7,Bpl7,Cpl7,Dpl7) ; 

[Aool8,Bool8, Cool8,Dool8] =series (Akouthigh, Bkouthigh, Ckouthigh,  Dkouthigh, 
Apl8,Bpl8,Cpl8,Dpl8) ; 

[Aool 9,  Bool  9 , Cool  9 , Dool 9 ] =series (Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh, 
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Apl9,Bpl9,Cpl9,Dpl9) ; 

[ Aoo2 0 , Boo2 0 , Coo2 0 , Doo2 0 ] =series ( Akouthigh, Bkouthigh, Ckouthigh,  Dkouthigh,  . 
Ap20,Bp20,Cp20,Dp20) ; 


%  Close  the  outer  loop  with  negative  feedback 

[ Aol , Bol , Col , Dol ] =cloop (Aool , Bool , Cool , Dool , -1 ) ; 

[ Ao2 , Bo2 , Co2 , Do2 ] =cloop ( Aoo2 , Boo2 , Coo2 , Dool , -1 ) ; 
[Ao3,Bo3, Co3,Do3]=cloop (Aoo3,Boo3,Coo3,Doo3,  -1)  ; 

[ Ao4 , Bo4 , Co4 , Do4 ] =cloop ( Aoo4 , Boo4 , Coo4 , Doo4 , -1 ) ; 
[Ao5,Bo5, Co5,Do5] =cloop (Aoo5,Boo5,Coo5,Doo5, -1) ; 
[Ao6,Bo6, Co6,Do6] =cloop (Aoo6, Boo6, Coo6, Doo6,  -1) ; 
[Ao7,Bo7, Co7,Do7] =cloop (Aoo7,Boo7,Coo7,Doo7, -1) ; 
[Ao8,Bo8, Co8,Do8] =cloop (Aoo8,Boo8,Coo8,Doo8, -1) ; 
[Ao9,Bo9, Co9,Do9] =cloop (Aoo9,Boo9,Coo9,Doo9,  -1)  ; 
[AolO,BolO, ColO,DolO] =cloop ( AoolO, Bool 0, Cool 0,  Dool 0, -1) ; 
[ Aol 1, Boll, Coll, Doll] =cloop ( Aooll, Bool 1, Cool 1,  Dool 1, -1) ; 
[Aol2,Bol2, Col2,Dol2] =cloop (Aool2,Bool2,Cool2,Dool2, -1) ; 
[Aol3,Bol3, Col3,Dol3]=cloop (Aool3,Bool3, Cool3,  Dool3,  -1) ; 
[Aol4,Bol4,Col4,Dol4]=cloop (Aool4,Bool4,Cool4,Dool4, -1) ; 
[Aol5,Bol5,Col5,Dol5]=cloop (Aool5,Bool5,Cool5,Dool5, -1) ; 
[Aol6,Bol6,Col6,Dol6]=cloop (Aool 6, Bool 6, Cool  6, Dool 6, -1) ; 
[Aol7,Bol7,Col7,Dol7]=cloop (Aool7,Bool7,Cool7,Dool7, -1)  ; 
[Aol8,Bol8,Col8,Dol8] =cloop (Aool8,Bool8,Cool8,Dool8, -1) ; 
[Aol9,Bol9,Col9,Dol9] =cloop (Aool9, Bool 9, Cool  9, Dool 9, -1) ; 
[Ao20,Bo20, Co20,Do20] =cloop (Aoo20,Boo20, Coo20,Doo20, -1) ; 

%  Get  step  response  of  closed  outer  loop 
t=0:.l:5; 

[yl,  xl] =step (Aol, Bol, Col, Dol, 1, t) ; 

[y2,  x2] =step (Ao2,Bo2, Co2,Do2, 1, t) ; 

[y3, x3] =step (Ao3, Bo3, Co3,Do3, 1, t) ; 

[y4,x4] =step (Ao4,Bo4, Co4,Do4,  l,t)  ; 

[y5, x5] =step (Ao5,Bo5, Co5,Do5, 1, t) ; 

[y6, x6] =step (Ao6,Bo6, Co6,Do6, 1, t) ; 

[y7, x7] =step (Ao7,Bo7, Co7,Do7, 1, t) ; 

[y8,x8]=step (Ao8,Bo8, Co8,Do8, 1, t) ; 

[y9, x9] =step (Ao9,Bo9, Co9,Do9, 1, t) ; 

[yl0,xl0] =step (Aol0,Bol0, Col0,Dol0,  l,t)  ; 
[yll,xll]=step(Aoll, Boll, Coll, Doll, l,t) ; 

[yl2,xl2] =step (Aol2,Bol2, Col2,Dol2, l,t) ; 

[yl3, xl3] =step (Aol 3, Bol 3, Col3,Dol3, 1, t)  ; 


B-18 


[yl4,xl4]=step (Aol4,Bol4, Col4,Dol4, l,t) ; 

[yl5, xl5] =step (Aol5,Bol5, Col5,Dol5, 1, t) ; 

[yl6, xl6] =step (Aol6,Bol6, Col6,Dol6, 1, t) ; 

[yl7,xl7]=step(Aol7,Bol7,Col7,Dol7,l,t) ; 

[yl8,xl8]=step(Aol8,Bol8,Col8,Dol8,l,t) ; 

[yl9,xl9]=step(Aol9,Bol9,Col9,Dol9,l,t) ; 

[y20, x20] =step (Ao20,Bo20, Co20,Do20, 1/ t) ; 

e=sum(abs (sum( [  (yl-yl3) ,  (yl-yl3) , (y2-yl3) , (y3-yl3) , (y4-yl3) , (y5-yl3) , (y6-yl3) , (y7-; 
(y8-yl3) , (y9-yl3) , (yl0-yl3) , (yll-yl3) , (yl2-yl3) , (yl4-yl3) , (yl5-yl3) ,  . . . 

(yl6-yl3) , (yl7-yl3) , (yl8-yl3) , (yl9-yl3) , (y20-yl3) ] ) ) ) 
figure  (6) 

plot (t,yl,t,y2,t,y3,t,y4,t,y5,t,y6,t,y7,t,y8,t,y9,t,yl0,t,yll,t,yl2, . . . 

t,yl3, t,yl4,t,yl5,t,yl6,t,yl7,t,yl8,t,yl9,t,y20) 

xlabel('Time  (sec)') 

ylabel ( 'Amplitude' ) 

grid 

axis{[0  5  0  1.3]) 
print  -deps  timebase 

%*******************************************  end  of  file  **** 

%*****************************************  file  relerr.m  **** 
function  j=relerr(]<) 

%  This  function  will  talce  the  given  minimal  order  Hinf 
%  controller  Keq  and  calculate  the  relative  error  to  the 
%  chosen  central  controller 

clear  numtemp  dentemp  nump  denp 
SV=[]  ;  j=0; 

F  =  -40;  Kf  =  [1  1] ;  G  =  .0247; 

1=1; 

for  i=l : 4 
Akeq(i)=F-G*k(l)  ; 

Bkeq(i, :)=Kf-G*[k(l+l)  k(l+2)]; 

Ckeq(i) =-k (1)  ; 

Dkeq(i, :)=-[k(l+l)  k(l+2)]; 

1=1+3; 

end 

Akcen=Akeq(3) ; 

Bkcen=Bkeq ( 3 , : ) ; 
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Ckcen=Ckeq(3) ; 

Dkcen=Dkeq(3, : ) ; 

%  Input  the  dynamics  at  the  20  flight  conditions 
numSYS=20; 

Along=[47.4,0; [-.2296  . 9931; . 02436  2046] ;  %Am3h26 

68.5,0; [-.2423  . 9964; -2 . 342  -. 1737] ;  %Am5h40 
100.1,0; [-.4285  . 9916; -. 7473  -. 3123] ;  %Am4h22  *3 

158.4,0; [-.5088  .994;-!. 131  -.2804];  %Am6h30 
189.9,0; [-.8018  .9847;-!. 521  -.5944];  %Am4h6 
255.0,0; [-.8930  . 9852;-4 . 1582  -. 6873] ;  %Am5hl0  *6 
301.1,0; [-.9181  . 9872;-6. 2419  -.6920] ;  %Am6hl5  *7 
355.0,0; [-.9920  .  9888;-7 . 8450  -. 7525] ;  %Am7hl8_5  *8 
426.4,0; [-1.175  .9871;-8.458  -.8776];  %Am7hl4 
496.0,0; [-1.4710  . 9808; -11 . 5022  -1 . 0846] ;  %Am6h2  *10 
557.0,0; [-2.1163  .  9872; -32 . 6459  -1 . 1826] ;  %Am8hl4  *11 
603.0,0; [-1.562  .9862;-14.94  -1.132];  %Am8hl2; 

614.4,0;  [-1.905  .9895;-33.88  -.9872];  %M95h20; 

652.0,0; [-1.675  .  9853;-16. 16  -1.212] ;  %Am8hlO; 

705.0,0; [-2.1163  .  9872; -32 . 6459  -1 . 1826] ;  %Am9hl4  *15 
789.1,0; [-1.994  .9828;-19.44  -1.427];  %Am8h5; 

825.2,0; [-2.452  . 9856;-38 . 61  -1 . 34] ;  %Am9hlO; 

890.8,0; [-2.328  .9831;-30.44  -1.493]  %Am85h5; 

956.0,  0; [-2.8375  .  9855; -51 . 8325  -1 . 4037] ;  %Am95h9  *19 

998.7,0; [-2.911  .9835;-46.47  -1.553];  %Am9h5; 
inf, 0]  ; 

C=[l  0;0  1]; 

D=[0;0]; 

Blong=[0;l]  ; 

w=logspace (-2,2,50); 
int=[250  500  750  1000] ; 
r=l ; i=l ; count=l ; 
for  i  =  1 : 4 

while  Along (r, 1)  <  int(i) 

[Ap,Bp, Cp,Dp] =feedback (Along ( r+1 :r+2, : ) ,Blong,  C,D, Akeq(i) , . . 
Bkeq(i, :) ,Ckeq(i) ,Dkeq(i, :) ,1) ; 

[nump(count, :) ,denp(count, :) ] =ss2tf (Ap,Bp, Cp (1, :) ,Dp(l, :) ) ; 
[svcla] =sigma(Ap,Bp,Cp (1, : ) ,Dp (1, : ) , w) ; 

SVcla= [SVcla; svcla] ; 

[svclp]=sigma(Ap,Bp,Cp(2, :) ,Dp(2, ;) ,w) ; 
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SVclp= [SVclp; svclp] ; 
r=r+3; 

count=count+l; 

end; 

i=i+l; 

end; 


%  Reduce  the  transfer  function  to  essential  parts 

%  (chop  off  added  zeros) 

nump=nump ( : , 3 : 4 ) ; 

numpO=nump (13,:); 

denpO=denp (13, :) ; 

%  Compute  Relative  Error  (Po-P)Po"-l 
for  z=l:20 

numtemp=conv (nump ( z , : ) , denpO) -conv (numpO , denp ( z , : ) ) ; 
dentemp=conv (denp ( z , : ) , numpO) ; 

[a,b, c,d] =tf2ss (numtemp,dentemp) ; 

[sv]=sigma(a,b,c,d,w) ; 

SV=[SV;sv] ; 
j=j+max (sv) ; 
end 

return 

%*******************************************  end  of  file  **** 

%*****************************************  file  relend. m  **** 
function  j=relend(k) 

%  This  function  will  talce  the  given  minimal  order  Hinf 
%  controller  Keq  and  calculate  the  relative  error  to  the 
%  chosen  central  controller 

numtemp=  [  ] ;  dentemp=  [  ] ;  niamSYS=0 ;  Ap=  [  ] ;  Bp=  [  ] ;  Cp=  [  ] ;  Dp=  [  ] ;  j=  [  0  ] ; 


n=45; 

F  =  -40;  Kf  =  [1  1] ;  G  =  .0247; 

1=4; 

for  i=l : 4 
Akeq(i) =F-G*k (1)  ; 
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Bkeq(i, :)=Kf-G*[k(l+l)  k(l+2)]; 

Ckeq(i) =-k (1) ; 

Dkeq(i, :)=-[k(l+l)  k(l+2)]; 

1=1+3; 

end 

%  Following  are  the  flight  conditions  used  in  the 
%  evaluation  of  relative  error. 

A=[47.4,0;  [-.2296  .  9931;  .02436  -.2046] ;  %Aiti3h26 
68.5,0; [-.2423  . 9964;-2. 342  -.1737] ;  %Am5h40 
100.1,0; [-.4285  .  9916; -. 7473  -. 3123] ;  %Am4h22 
158.4,0; [-.5088  .994;-!. 131  -.2804];  %Am6h30 
189.9,0; [-.8018  . 9847; -1 . 521  -. 5944] ;  %Am4h6 
255.0,0; [-.8930  . 9852; -4 . 1582  -. 6873] ;  %Am5hl0 
301.1,0; [-.9181  .9872;-6. 2419  -.6920];  %Am6hl5 
355.0,  0;  [-.9920  .  9888;-7. 8450  -.7525] ;  %Ain7hl8_5 
426.4,0; [-1.175  .9871;-8.458  -.8776];  %Am7hl4 
496.0,  0;  [-1.4710  .  9808;-ll  .5022  -1.0846];  %M6h2 
557.0,  0; [-2.1163  .  9872;-32 . 6459  -1.1826];  %Am8hl4 
603.0,0; [-1.562  .9862;-14.94  -1.132];  %Am8hl2; 
614.4,0; [-1.905  . 9895; -33 . 88  - . 9872] ;  %Am95h20; 
652.0,  0;  [-1.675  .  9853; -16 . 16  -1 .212] ;  %Ain8hlO; 
705.0,  0;  [-2.1163  .  9872;-32 . 6459  -1.1826];  %Ain9hl4 
789.1,0; [-1.994  .9828;-19.44  -1.427];  %Am8h5; 
825.2,0; [-2.452  . 9856; -38 . 61  -1 .34] ;  %Am9hlO; 
890.8,0; [-2.328  .9831;-30.44  -1.493]  %Am85h5; 
956.0,0; [-2.8375  . 9855;-51 . 8325  -1.4037];  %Am95h9 
998.7,0; [-2.911  .9835;-46.47  -1.553];  %Am9h5; 
inf, 0]  ; 

C=[l  0;0  1]; 

D=[0;0]; 

Blong= [0; 1]  ; 

%  note  i  =  interval  number 

int=[k(l:3)  1000]; 
s=length (int) ; 

for  1=1: (s-1) 
for  i=  1 : (s-1) 
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if  int (i) >int (i+1) 
temp  =  int (i) ; 
int (i)  =  int (i+1) ; 
int (i+1) =temp; 
end 
end 
end 


r=l;i=l; 
while  i  <=  4 
while  A(r, 1)<  int(i) 
numS  Y  S =numS  Y  S + 1 ; 

[Ap,Bp, Cp,Dp] =feedback (A(r+1 : r+2, : )  ,Blong, C,D,  . . . 

Akeq(i) ,Bkeq(i, : ) , Ckeq(i) ,Dkeq(i, : ) ,  1) ; 

[nump (numSYS, : ) , denp (numSYS, : ) ] =ss2tf (Ap,Bp, Cp (1, : ) ,Dp (1, : ) ) ; 
r=r+3; 
end; 
i=i+l; 
end; 

w=logspace(-2,2,n) ; 

numpO=nump (13,  : )  ; 
denpO=denp (13,  : )  ; 

for  z=l: numSYS 

numtemp=conv(nump  (z,  : )  ,  denpO)  -conv(niampO,  denp  (z,  : ) )  ; 
dent emp=conv (denp (z, : ) , numpO) ; 
numtemp=numtemp ( : , 3 : 7 ) ; 
dentemp=dentemp ( : ,  3 : 7 ) ; 

[a,b, c, d] =tf2ss (numtemp,dentemp) ; 

[sv]=sigma(a,b,c,d,w)  ; 
j= j+max (sv) ; 
end 

return 

%*******************************************  end  of  file  **** 

%*****************************************  file  relnum.m  **** 
function  j=relnum(k) 
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%  This  function  will  take  the  given  minimal  order  Hinf 
%  controller  Keq  and  calculate  the  relative  error  to 
%  the  chosen  central  controller 

numtemp= [ ] ; dentemp= [ ] ; numSYS=0 ; j= [ 0 ] ; 

n=45; 

F  =  -40;  Kf  =  [1  1] ;  G  =  .0247; 

%  Calculate  the  controller  for  each  interval 
%  Note  i  =  interval  number 


1=10; 

for  i=l :k  (1) 

Akeq(i)=F-G*k(l) ; 

Bkeq(i, :)=Kf-G*[k(l+l)  k(l+2)]; 

Ckeq(i)=-k(l)  ; 

Dkeq(i, :)=-[k(l+l)  k(l+2) ] ; 

1=1+3; 

end 

%  Following  are  the  flight  conditions  used  in  the  evaluation 
%  of  relative  error, 

A=[47. 4,0; [-,2296  . 9931; . 02436  -.2046] ;  %Am3h26 

68.5,0; [-.2423  , 9964; -2 . 342  - . 1737] ;  %Am5h40 
100,1,0; [-,4285  . 9916; -. 7473  -. 3123] ;  %Am4h22 
158.4,0; [-,5088  . 994;-l .131  -.2804] ;  %Am6h30 
189.9,0; [-.8018  .9847;-!, 521  -.5944];  %Am4h6 
255.0,0; [-.8930  . 9852; -4 . 1582  -. 6873] ;  %Am5hl0 
301.1,0; [-.9181  . 9872;-6. 2419  -.6920] ;  %Am6hl5 
355.0,0; [-.9920  . 9888; -7 , 8450  - .7525] ;  %Am7hl8_5 
426.4,0; [-1.175  .9871;-8.458  -.8776];  %Am7hl4 
496.0,0; [-1.4710  . 9808;-11.5022  -1.0846];  %Am6h2 
557.0,0; [-2.1163  .  9872;-32 . 6459  -1.1826];  %Am8hl4 
603.0,0; [-1.562  .9862;-14.94  -1.132];  %Am8hl2; 

614.4,0; [-1.905  .9895;-33.88  -.9872];  %Am95h20; 

652.0,0; [-1,675  . 9853; -16 . 16  -1 .212] ;  %Am8hlO; 

705.0,0; [-2.1163  . 9872; -32 . 6459  -1.1826];  %Am9hl4 
789.1,0; [-1,994  ,9828;-19.44  -1.427];  %Am8h5; 

825.2,0; [-2.452  , 9856; -38 . 61  -1 . 34] ;  %Am9hlO; 

890.8,0; [-2.328  .9831;-30.44  -1.493]  %Am85h5; 
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956.0,  0; [-2.8375  .  9855; -51 . 8325  -1.4037];  %Am95h9 
998.7,0; [-2.911  .9835;-46.47  -1.553];  %Mi9h5; 
inf, 0] ; 

C=[l  0;0  1]; 

D=[0;0]  ; 

Blong=[0; 1]  ; 

%  Sort  the  interval  break  points  from  lowest  to  highest 

int=[k(2:k(l) ) ,1000] ; 
s=length (int) ; 

for  1=1: (s-1) 
for  i=  1 : (s-1) 
if  int (i) >int (i+1) 
temp  =  int (i) ; 
int (i)  =  int (i+1) ; 
int (i+1) =temp; 
end 
end 
end 

%  Check  the  dynamic  pressure  of  the  flight  condition  and 
%  close  the  inner  loop  with  the  appropriate  controller 
%  for  the  interval. 

%  Then  convert  from  state  space  to  transfer  function  form. 

r=l;i=l; 
while  i  <=  k (1) 
while  A(r,l)<  int(i) 
numS Y  S =numS  Y  S + 1 ; 

[Ap,Bp,Cp,Dp]=feedback(A(r+l;r+2, :) ,Blong,C,D,  . . 

Akeq(i) ,Bkeq(i, : ) ,Ckeq(i) ,Dkeq(i,  : )  ,  1)  ; 

[nump (numSYS, : ) , denp (numSYS, :) ] =ss2tf (Ap,Bp, Cp (1, :) ,Dp(l, :) ) ; 
r=r+3; 
end; 
i=i+l; 
end; 

%  Select  the  central  controller  and  calculate  the  relative  error 
w=logspace (-2,2,n) ; 
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numpO=nump (k  (37) ,  : ) ; 
denpO=denp (k (37) , : ) ; 

for  z=l:nuinSYS 

numt emp=conv ( nump ( z , : ) , denp 0 ) -conv ( numpO , denp ( z ,  : ) )  ; 
dentemp=conv ( denp ( z , : ) , numpO) ; 
numtemp=nunitemp  ( : ,  3 : 7 )  ; 
dentemp=dentemp ( : , 3 : 7 ) ; 

[a,b, c, d] =tf2ss ( numt emp, dent emp) ; 

[sv]=sigma(a,b,c,d,w) ; 
j= j+max (sv) ; 
end 

j=j 

return 

%*******************************************  end  of  file  **** 

%***************************************  file  relinter.m  **** 
function  j=relinter (k) 

%  This  function  will  take  the  given  minimal  order  Hinf 
%  controller  Keq  and  calculate  the  relative  error  to 
%  the  chosen  central  controller 


e= [ ] ; numSYS=0 ; j= [ 0 ] ; 


n=45; 

F  =  -40;  Kf  =  [1  1] ;  G  =  .0247; 

%  Following  are  the  flight  conditions  used  in  the 
%  evaluation  of  relative  error. 

A=[47.4,0; [-.2296  . 9931; . 02436  -.2046] ;  %Am3h26 

68.5,0; [-.2423  . 9964; -2 . 342  - . 1737] ;  %Am5h40 
100.1,0; [-.4285  .  9916;-. 7473  -.3123] ;  %Am4h22 
158.4,0; [-.5088  . 994; -1 . 131  - .2804] ;  %Am6h30 
189.9,0; [-.8018  . 9847;-! . 521  -. 5944] ;  %Am4h6 
255.0,0; [-.8930  . 9852;-4 . 1582  -. 6873] ;  %Am5hl0 
301.1,0; [-.9181  . 9872;-6. 2419  -.6920] ;  %Am6hl5 
355.0,0; [-.9920  . 9888;-7 . 8450  -.7525] ;  %Am7hl8_5 
426.4,0; [-1.175  . 9871; -8 . 458  - . 8776] ;  %Am7hl4 
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496.0,0; [-1.4710  . 9808;-ll .5022  -1.0846];  %Am6h2 
557.0,0; [-2.1163  .  9872;-32 . 6459  -1.1826];  %Am8hl4 
603.0,0; [-1.562  .9862;-14.94  -1.132];  %Am8hl2; 

614.4,0; [-1.905  . 9895;-33 . 88  -. 9872] ;  %Ain95h20; 

652.0,0;  [-1.675  .  9853;-16. 16  -1.212] ;  %Ain8hlO; 

705.0,0; [-2.1163  . 9872;-32 . 6459  -1.1826];  %Am9hl4 
789.1,0;  [-1.994  .9828;-19.44  -1.427];  %Ain8h5; 

825.2,0; [-2.452  . 9856; -38 . 61  -1 . 34] ;  %Am9hlO; 

890.8,0; [-2.328  .9831;-30.44  -1.493]  %Am85h5; 

956.0,0; [-2.8375  . 9855; -51 . 8325  -1.4037];  %Am95h9 
998.7,0; [-2.911  .9835;-46.47  -1.553];  %Am9h5; 
inf, 0] ; 

C=[l  0;0  1]; 

D=[0;0]; 

Blong=[0; 1] ; 

%  Sort  the  interval  break  points  from  lowest  to  highest 

int=[45,k(2:k(l) ) ,1000] ; 
s=length (int) ; 

for  1=1: (s-1) 
for  i=  1 : (s-1) 
if  int (i) >int (i+1) 
temp  =  int (i) ; 
int  (i)  =  int (i+1)  ; 
int (i+1) =temp; 
end 
end 
end 

%  Check  the  dynamic  pressure  of  the  flight  condition  to 
%  determine  which  interval  endpoints  to  use.  Calculate 
%  the  corresponding  values  of  N  and  M. 

%  Form  the  controller  Keq,  close  the  inner  loop  and 
%  convert  to  transfer  function  form  for  the  relative  error 
%  calculation. 

ni=10;mli=20;m2i=30; 
r=l; i=l; 
while  i  <=  k (1) 


while  A(r,l)<  int(i+l) 
numSYS=numSYS+l ; 

N=k (ni) + (k (ni+1) -k (ni) ) * ( (A(r, 1) -int (i) ) / (int (i+1) -int (i) ) ) ; 

Ml=k (mli)  +  (k (mli+1) -k (mli) ) * ( (A(r, 1) -int (i) ) / (int (i+1) -int (i) ) )  ; 
M2=k  (m2i)  +  (k  (in2i+l)  -k  (m2i) )  *  ( (A(r,  1)  -int  (i) )  /  (int  (i+1)  -int  (i) ) )  ; 
Akeq=F-G*N; 

Bkeq=Kf-G* [Ml  M2]; 

Ckeq=-N; 

Dkeq=-[M1  M2]  ; 

[Ap,Bp, Cp,Dp] =feedback (A(r+1 : r+2, : ) ,Blong, C,D,  . . . 

Akeq, Bkeq, Ckeq, Dkeq, 1 ) ; 

[nump (numSYS, : ) , denp (numSYS, : ) ] =ss2tf (Ap,Bp, Cp (1, : ) ,Dp (1,  : ) )  ; 
r=r+3; 
end; 

i=i+l;ni=ni+l;mli=mli+l;m2i=m2i+l; 

end; 

%  Select  the  central  controller  and  calculate  the  relative  error 

w=logspace(-2,2,n) ; 
nump=nump ( : , 3 : 4 ) ; 
numpO=nump (k (40) , : ) ; 
denpO=denp (k (40) , : ) ; 

for  z=l; numSYS 

numt  emp=conv ( nump ( z , : ) , denp  0 ) -conv ( nump  0 , denp ( z , : ) ) ; 
dentemp=conv ( denp ( z , : ) , numpO) ; 

[a,b, c, d] =tf2ss (numtemp, dentemp) ; 

[sv]=sigma(a,b,c,d,w) ; 
j= j+max (sv) ; 
end 

%  Close  the  outer  loop  and  calculate  the  step  response 
longout  %  a  file  containing  the  outer  loop  controller 
for  z=l:5 

[nums(z, :) ,dens(z, :) ] =series (numoutlow, denoutlow,  . . . 
nump ( z , : ) , denp ( z , : ) ) ; 

[numcl  (z,  : ) ,  dencl  (z,  ; )  ]  =cloop  (niams  (z, : ) ,  dens  (z,  : ) ,  -1)  ; 
end 

for  z=6: numSYS 
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[nums (z, : ) , dens (z, : ) ] =series (numouthigh, denouthigh, . . . 
nump(z,  :)  ,denp(z,  :))  ; 

[numcl(z, :) ,dencl(z, :) ] =cloop (nums (z, :) ,dens(z, :) ,-l) ; 
end 

t=0: .1:5; 


for  z=l:numSYS 

[y, x] =step (numcl (z, : ) , dencl (z, : ) , t) ; 

Y=[Y,y];X=[X,x]; 

end 

%  Calculate  the  error  of  the  step  response 

for  z=l:numSYS 
e(:,z)=Y(:,z)-Y(:,k(40)); 
end 

e=sum(abs (sum(e) ) )  ; 

j=j+e 

return 

%*******************************************  end  of  file  **** 

%****************************************  file  longout.m  **** 

%  Outer  loop  longitudinal  controller 

Kouthigh=l . 1995e2; 

zouthigh= [-2 . 418e3, -5 .4039+7.25721,-5 .4039-7.25721]'; 

pouthigh= [-7 . 1836el+4 . 0476eli,-7 . 1836el-4 . 0476eli, -7 .2319e-2, -2 . 0509el] ' ; 
Koutlow=l .2951e2; 

zoutlow=[-6.0168e2, -3. 8407+6. 25071, -3. 8407-6. 25071] ' ; 
poutlow=[-5.9135el+5.1963eli,-5.9135el-5.1963eli,-3.6321e-2,-6.6426]' ; 

[numouthigh,  denouthigh] =zp2tf (zouthigh,pouthigh,Kouthigh) ; 

[numoutlow, denoutlow] =zp2tf (zoutlow,poutlow,Koutlow) ; 

[Akoutlow, Bkoutlow, Ckoutlow, Dkoutlow] =zp2ss ( zoutlow, poutlow, Koutlow) ; 

[  Akouthigh, Bkouthigh, Ckouthigh, Dkouthigh] =zp2ss ( zouthigh, pouthigh, Kouthigh) ; 

%*******************************************  end  of  file  **** 
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